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1. DIGITAL APPROXIMATION BY POINT BY POINT STATE COMPARISON 
1.1 Introduction 

The problem is that of approximating a continuous-data system 
by a sampled-data model by comparison of the states of the two 
systems. 

Consider that the continuous-data system as shown in Fig. l - i 
is described by the following time- invariant dynamic equations: 


x^(t) = Ax^(t) + Bu(t) (1-1) 

u(t) = E(0)r(t) - GCOjiMt) (1-2) 

where 

x (t) = n x 1 state vector 
— c 

£(t) = m x 1 control vector 
r(t) = m x 1 input vector 
A = n x n coefficient matrix 
B = n x m coefficient matrix 
E (0) = I = m x m identity matrix 
G(0) = m x n feedback matrix 
The initial state is given by ^(0) = 

Substituting Eq. (1-2) into Eq. (1-1) yields 

x (t) = [A - BG (0) 1 x (t) + BE (0) r (t) (1-3) 

— C —C~ 

The solution of Eq. (1-3) for t > t is 

— o 


x_(t) =■ 


[A-BG(O) ] (t-t.) 
e Ox 


<V 


♦r 


(A-BG (0) ] (t-T) 


BE (0) r (T) dt (1-4) 


o 
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where 

[A-BG(O) ] (t-t ) “ • 

e u = y 4- [A - BG(0)] J (t - t) J (1-5) 

, L n n 1 0 

3=0 J 

The block diagram of the sampled-data system which is to 
approximate the system of Fig. 1-1 is shown in Fig. 1-2. The outputs of 
the sample-and-hold devices are a series of step functions with 
amplitudes denoted by u^kT) for kT <_ t < (k + 1)T. The notation, 

G(T) and E(T) , denote the feedback gain and the forward gain of the 
sampled-data system, respectively. The dynamic equations for the 
sampled-data model are : 

x (t) = Ax (t) + Bu (kT) x (0) = x (1-6) 

— s — s — s — s “SO 

u (kT) = E(T)r(kT) - G(T)x (kT) (1-7) 

— s — — s • 

for kT £ t £ (k + 1)T. 

The A and B matrices in Eq. (1-6) are identical to those of Eq. (1-1) 
Substituting Eq. (1-7) into Eq. (1-6) yields 


x (t) = Ax (t) + B[E(T)r(kT) - G(T)x (kT) ] 


( 1 - 8 ) 


for kT £ t £ (k + 1)T. 

The solution of Eq. (1-8) with t = (k + 1)T and t = kT is 

o 


. (v+i ) t 

r/, I at f A(kT+T-T) . 

x_ [ (k + 1)T] = | e - | e dxBG(T) 

kT 


x (kT) 
— s 


r(k+l)T 


kT 


A(kT+T-T) . . 
e dx,BE(T)r (kT) 


(1-9) 
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The problem is to find E (T) and G(T) so that the states of the 
sampled-data model are as close as possible to that of the continuous - 
data system at the sampling instants, for a given input r^(t) . Further- 
more, in order that the solution for E(T) is independent of r^(t) it 
is necessary to assume that £(T) => r^(kT) for kT £ t<(k + 1)T. 

Therefore, effectively, the input of the continuous-data system of 
Fig. 1-1 is assumed to pass through a sample-and-hold device. The 
above assumption would not affect the solution if r^(T) has step 
functions as its elements. However, if the inputs are other than 
step functions , the approximation is a good one for small sampling 
periods. 

Now letting t = kT and t = (k + 1)T in Eq. (1-4) , and 
assuming r(x) = .r(kT) over one sampling period, we have 


x [(k + 1) T] 
— c 


[A-BG(O) ]T 
e 


x (kT) + 


(k+l)T 

kT 


[A-BG(O) ] (kT+T-T) 
e 


dTBE(0)r(kT) 


( 1 - 10 ) 


kT < t < (k + 1) T. 


The responses of Eq. (1-9) and Eq. (1-10) will match at 
t = (k + 1)T for an arbitrary initial state x^_(kT) and an arbitrary 
input r_(x) , if and only if the following two equations are satisfied. 


[A-BG(O) ]T AT 
e = e 


(k+1) T A[kT+T _ T] 
e 

kT 


dXBG(T) 


( 1 - 11 ) 




and 
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First working with Eq. (1-11), we let X = (k + 1)T - X. Then, 
Eq. (1-11) becomes 

e [A-BG(O) ]T = /T _ g AX ^ BG(T) (1 _ 13) 

Jo 

In principle, the feedback matrix, G (T) , of the sampled-data system 
can be determined from Eq. (1-13). However, it is subject to the 
limitations as discussed in the following section. 
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l . 2 Closed Form Solutions for G(T) 

In order that all n states of the continuous system x^OtT) match 
those of the sampled system x^CkT) at the end of each sampling 
period it is sufficient that equation (1-13) be satisfied. Equation 
(1-13) may be written in the simplified matrix form, 

D(T) = - 0(T)G(T) (1-14) 

where 

(1-15) 


9<t) 


f 1 fll 

e BdA 

0 


(1-16) 


Equation (1-14) consists of n scalar equations with mn unknowns. 
If the number of unknowns equals the number of equations, m = n, and 
0(T) is nonsingular, then a unique solution of Eq.(l-14) exists and 
is given by 


G(T) 


0 1 (T) D (T) 


(1-17) 


A Special Case for Exact Matching 

If the control system has more states than controls, n > m, which 
is the case in most control systems, Eq.(l-14) will generally not have 
a solution. However, if the rank conditions described in the 
following are satisfied, the system of equations in Eq.(l-14) are 
consistent and there still is a solution. 

Let 


G(T) 

= C V ^2 9 n ] 

(1-18) 

D(T) 

= [d lf d 2 , ... , d n ] 

(1-19) 
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where g. , i = 1, n, are m- dimensional vectors and d^, i=l, • ••/ 

n, are n-dimensional vectors, then if 

rank [0] = rank [0, dj (1-20) 

for all i=l, . .., n, the system of equations in Eq.(l-14) has at 
least one solution. If the above conditions are not satisfied the 
equations are inconsistent and no solution exists. 

Partial Matching of States 

Generally speaking, the above rank conditions will not be 
satisfied and thus for the case where n > m not all of the states 
of the continuous and sampled systems can be made to match at the 
end of each sampling period. 

Although it is not possible to match all of the states it can 
be shown that it is possible to match some of the states or algebraic 
sums of the states at each sampling period. Multiplying both sides 
of Eq. (1-14) by a constant m x n matrix H gives, 

HD(T) = - H0(T)G(T) (1-21) 

The above equation consists of mn scalar equations and mn 
unknowns. If H is chosen such that H0 (T) is nonsingular, Eq. ( 1 - 21 ) 
may be solved for a solution, ^(T) , 

G (T) - - [H0(T) ] _1 HD(T) (1-22) 

w 

It should be noted that the above solution for G (T) does not satisfy 

Eq. (1-14) . The reason is that when Eq. (1-14) is premultiplied by the 

2 

matrix H, it reduces the system of n equations to a system of mn 
equations. Only if Eq.(l-14) had been multiplied by a nonsingular 
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matrix would the new system of equations be equivalent to the 
original system. Since H is a m x n matrix, for n > m the 
transformation is a singular one, and the solution may not necessarily 
satisfy the original set of equations. 

To understand the consequencies of the transformation and the 
physical meaning of the solution of Eq.(l - 22), we match Eqs. (1-9) 
and (1-10) in the following form: 

x [ (k + 1) T] = <f> (T) x (kT) + 0 (T)E(0)r(kT) 

— C C "C C — 


where 


= x [(k + 1) T] = [<J> (T) " 0 (T)G(T) ]x. (kT) + 0 (T) E (t) r (kT) (1-23) 
— s 0 s 


. [A-BG (0) ]T 

(f> (T) = e 


(state transition matrix of the 
closed- loop continuous -data system) (1-24) 


AT 

4> 0 ( T) = e 


(state transition matrix of the 
open- loop system) 


0 (T) 
c 


• T e [ A -BG<0)]X aAB 
'0 


(1-25) 

(1-26) 


Premultiplying both sides of Eq. (1-23) by the m x n matrix H gives 


Hx [(k + 1) T] = H(j) (T) x (kT) + H0 (T)E(0)r(kT) 

— c c — c c 

= Hx [(k + 1) T] = H[cf> (T) - 0(T)G(T)]xJkT) + H0 (T) E (T) r (kT) 

S 0 is 

(1-27) 

For arbitrary x (kT) , x (kT) , and r (kT) , the last equation leads to 
c s 

H4» c (T) x c (kT) = H[4» 0 (T) - e(T)G(T)x g (kT) 


(1-28) 
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and 

H0 (T)E(0)r(kT) = H0(T)E(T)r (kT) (1-29) 

c — 

The significant point is that the solution, G^fT) , of Eq. (1-22) 

satisfies Eq- (1-29) for any arbitrary "initial" state x^ (kT) . 

Premultiplying both sides of Eqs- (1-9) and (1-10) by H simply 

transforms the n-dimensional state vectors x [ (k + 1)T] and 

— c 

x^[(k + 1)T] into a new m-dimensional vector, yj (k + 1)T] , such that 

y [ (k + 1) T] = Hx [k + 1) T] = Hx [ (k + 1)T] (1-30) 

c s 

Note that the new states y i [ (k + 1)T] , i = 1, 2, m, are algebraic 

sums of the original n state variables; that is, 


+ 


n 


1) T] = I 
j=l 


h . .x , [ (k 
13 3 


+ 1) T] 


i - I# 2, 


-m 


(1-31) 


Although the solution, G w (T), of Eq. (1-22) does not match the n 

states of the continuous system, x [ (k + 1)T] , to the n states of 

— c 

the sampled system x^[(k + 1)T] > at time (k + 1)T, it does match an 
algebraic sum of these states. 


£[(k + 1) T] = Hx [(k + 1 j T] = Hx [(k + 1) T] 


(1-32) 
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1.3 An Exact Solution for G (T) 

W 

With the understanding of the conditions and limitations as 

described in the preceding section, we shall derive an exact closed- 

form solution of the feedback gain G (T) which partially matches the 

w 

states at the end of each sampling period as governed by the weighting 
matrix H. Starting with Eq. (1-13) 


I 


Let 


then 


[ A-BG ( 0 ) ] T 
e 


AT 

e 


e A * dXBG(T) 


e«> (j + 


«, A j T j+1 


B 


[A-BG (0) ]T 
e 


= e AT - 0 (T) G (T) 


(1-33) 


(1-34) 


(1-35) 


Premultiplying by a weighting matrix H which is chosen such that 
H0(T) is nonsingular and solving for G(T) 


G (T) - [H0(T)3 _1 H[e AT - e t A “ BG(0)]T ] 
w 


(1-36) 


A special case of interest is when H = B ‘ . In this case H0 may 
be expressed as a quadratic form 


f T 


H0(T) = 


B'e A *B dX 


(1-37) 


Since 


- 00 3 3 

AA v A J T J 

e = i 

j=0 


jl 


(1-38) 
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AA 

for X sufficiently small e >0 (i.e. a positive definite matrix). 
Thus it can be argued that there exists a number a > 0 such that if 
T < a then 

e^> 0 for 0 < X < T 
and therefore the quadratic form 


t T 


aX 

B'e B dX > 0. 


Then for T < Ct, H0(T) is nonsingular. 
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1.4 G(T) By Series Expansion 

In this section we shall derive a Taylor series expansion for 
G(T) about T = 0. In general, if this series converges, G(T) can 
be approximated by talking a finite number of terms of the series 
expans ion. 

Let G(T) be represented by a Taylor series about T = 0; that is 


where 


G(T) 


x i t ' \ -l 

= lim G (T) = lim J — G 1 ^ (T)T 3 
K->oo Kt* 30 j=0 3 


(1-39) 


(T) = 9^G(T) 


(1-40) 


3t j 


T=0 


Substitute G (T) from Eq. (1-39) into Eq. (1-13) , we get 


[A-BG (0) ] T AT 
e = e 


f T 


aA . 


e dXB £ — G (D (T)T 3 (1-41) 

0 j=0 3 


Or 


X 

j=0 


[A - BG(0) l^T 3 

j! 


” M +1 


oo G, (T)T 
k 


V AT V AT Y k 

= jio + k £ 0 kI 


(1 


The last equation is written 


^ [A - BG(0) l^T 3 A^T* 1 A ] B 

v j ! j I < j + 1) ! 


X 

j=0 


00 _(k) /m ._k+j+l 

y G_ (T) T 


k=0 


k! 


j -0 


Now equating the coefficients of T^ (p = 1, 2, ...) to zero, we have 

[A - BG(0)] P _ A^ + P V 1 A P ~3~ 1 BG (3) (T ) = Q 

pi Pi jf 0 (p-j)ljl ~ U--44) 


-42) 


(1-43) 
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In general, it is possible to express G ^ P (T) in terms of ^ (T) , 

G (P " 3) (T) , , G (1) (T) , G (0) (T) , where G (0) (T) = G(0) . 


Equation (1-44) is written 


[A-BG(O) ] P _ A^ P r 2 A^^BG^ (T) BG (p 1) (T) 
pi " pi + _.“ 0 (p - 1) l j I (p - 1) I 


(1-45) 


Therefore, 


BG^ P ~ X ^ (T) 
(p - 1) I 


A_ 

P! 


[A-BG(0)1 P Y A m ~^ 1 BS (3) (T) 




j=o 


(P - j) ! jl 


(1-46) 


Then, let H be an m x n matrix such that (HB) is nonsingular; Eq. (1-46) 


(P — 1) 

leads to the solution of G (T) , 


G (P 1} (T) « (HB) -1 H 


Ar 

(P 


[A - BG (0) ] J 
P 


p — 1 f 2/ • • • 


P-2 P~j-1 R (j) , J 

(p - 1) J l P5 .-- LT) (1-47) 

* j _ 0 <P “ 3) 131 j 


(p-1) 


The following table gives the expressions for G (T) for 


p = 1, 2, and 3. Then, according to Eq. (1-39), 


G (T) = G (0) (T) = G (0) 

(1-48) 

G £ (T) = G(0) + TG (1) (T) 

(1-49) 

111 T 2 (2) 

G 3 (T) = G(0) + TG (T) + — G (T) 

(1-50) 

Table 1-1. 



P 

G (P-1) (T). 

1 

G (0) 

2 

“ G (0) {A - BG(0) ] 

3 

(HB)“ H{- -- A3G(0) [A - BG(0) ] + j BG (0) [ A - BG (0) ] } 
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It is interesting to note that since (T) and G (T) are 

not dependent upon H, the one-term and the two-term series 
approximations of G(T) , G^T) and G 2 (T) , respectively, will attempt 
to match all the states of the continuous-data and the sampled-data 

! 

^systems. Beyond two terms , the weighting matrix, H, must be used, 

i 

! and only certain states and combinations of states are matched, 
depending on the H selected. As more terms are added to the series 
approximation, G K (T) will approach the exact solution of (MT) of 


Eq. (1-36). 



15 


1.5 An Exact Solution For E (T) 

The solution of E(T) is obtained from Eq. (1-23) < 

6 (T) E (T) = 6 (T) E (0) (1-51) 

c 

Similar to the closed-form solution for G(T) discussed in 
Section 1.2, if m = n, and 0(T) is nonsingular, a unique solution of 
Eq. (1-51) exists, and is given by 

E (T) = 6 _1 (T)e (T)E(O) (1-52) 

C 

When n > m, the solution of the matrix E(T) which corresponds 
to the partial matching of states is obtained from Eq. (1-29) , 

H0(T)E (T) = H0 (T)E(O) (1-53) 

w c 

Therefore , 

E (T) = [H0 (T) ] _ 1 H0 (T)E(O) (1-5A) 

w c 

where it is assumed that H0(T) is nonsingular. 


<> 
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1 . 6 E(T) By Series Expansion 

Similar to the solution of G(T) , the matrix E (T) may be 
expanded into a Taylor series about T = 0; that is 


E (T) =lim E (T) 
K-x» K 


where 


K-l ... 

E (T) = l 77 E°' <T)T 3 
j=0 3 


where 


E lj) (T) - 


3t- 


T=0 


(1-55) 


(1-56) 


(1-57) 


Substituting the series expansion of E (T) into Eq. (1-51), we get 


0 (T) l E (3) (T)T 3 = 0 (T)E(O) 

j=o 3 c 


(1-58) 


Or, 


oo i +1 00 

i J t 3 J:„ v i_„oo 


,1 r « uu <« k ■ { fr » - »«>»’ fn 


■i t 3+1 

I 3 4 BE (0) (1-59) 


The last equation is simplified to 


l 

5-0 


Ah. ? E (kl (TIT^-" 1 

<0 + 11 1 yto 


[A - 


BG(0) ] 3 T 3+1 


k! 


(j + 1) ! 


BE (0) 


) 


= 0 (1-60) 
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Equating the coefficients of (p = 1, 2, ...) to zero, we have 


Or, 


[A - BG(0) ] P_1 BE(0) . A P_1 B A P_2 BE (1) (T) 

+ - E(0) + 7— 7TE 

p! pi (p - 1) I 


^ p " 3 M . ( 2 . V t ... t aV'p- 1 ;^ . u . 61) 


(p - 2) 12! 


1! (p - 1) 1 


BE (P 1) (T) = 


[A - BG{0) ] P_1 BE(0) , P r 2 A P j 1 BE (j) (T) 

p (p " 1} 1 (P - "jTTJT 


(1-62) 


Now let H be an m x n matrix such that HB is nonsin gular , the solution 
of E vp ; (T) is given by 


E (p 1) (T) = (HB) -1 h{ 


[A - BG(0) ] P-1 BE(0) 


. ... P Y 2 A P - j - 1 BE (j) (T ) 1 n ... 

" (p " 1)1 1 I ; (1 ' 63) 


j=0 


p — 1, 2, ••• 

The following table gives the expressions for E^ P ^ (T) for 


p = 1 , 2 , and 3 . 


Table 1-2. 


p 

E (P-1) (T) 

1 

E (0) 

2 

- ^ G (0) BE (0) 

3 

(hbi-M “f 013 - “<0)»- BG(0)1B } E «» 

^6 3 




Notice that, similar to the situation for G (T) , if only up to 
two terms are used in the series approximation for E(T) , the H matrix 
is not needed, and all the states of the sampled-data and the continuous 
data systems are apparently matched. 
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1.7 Computer Program for Computation of Truncated Series, G ^ (T) and E ^JT) 
The power series expansions of G^(T) and E w (T) are defined in 
Eqs. (1-36) and (1-54), respectively. An approximation utilizing the 
first K terms of the infinite power series can be written as 


K-i ... T i 

G (T) = G (T) = l G V (T) — 
W K i=0 11 


(1-64) 


K-l m T i 

and E (T) = E (T) = £ E V '(T) — 

w * - i==0 


where 


(1-65) 


i - 2 (t) j 


G U ‘ 1) (THH»,' 1 K ^ j o (i . jMjl 

(A- 3) 

• i-2 , i — j — 1 , . . 

-l„r(A - BG(0)3' 1_1 BE(O) - (i-l)i l * i T BE 3 < t: 

"l i j=0 - J 


E (l 1) (T)=<HB) 


G (0) = G(0) 


(1-67) 

( 1 - 68 ) 


and E^ = E(0) 


(1-69) 


Equation (1-64)- (i- 67 )ca n be modified for computational efficiency, 
Let 

6 (m) ( T ) = ^ (1-70) 


(m - 1) I 


E (m - 1) 1 


(1-71) 


1 - 66 ) 

} 
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i,-4 


( 1-72 ) 


[A - BG ( 0) ] 


[A - BG (0) ] ± = ±J 


(1-73) 


(AB) 


A 1 B 


i (i + 1) 1 


(1-74) 


then the K-term series approximations of G (T) and E (T) become 

w w 


and 


G (T) = l G U+1) (T)T 1 
K i=0 


v >s(i+l) i 
E (T) = T £ 1+1 (T) T 

K i-0 


(1-75) 


(1-76) 


with 


»(i) 


- 1 . 


e vx/ (T) = (HB) H 


A. - [A-BGCOjrl (AB). ... G { ^ _1 ^ (T) 
i i . o. i-i+i 


j=2' i-3 + l 


( 1-77 ^ 


(t) = (HB) *" 1 H 


(A-BG(Q) ] i _ 1 BE(0) 


3=2 


e' 3 ' 1 ’ (T) 


(1-78) 


e (1) (T) = G (0) 


and £ (1) (T) = E(0) 


(1-79) 

(1-80) 


Figure 1-3 shows the flow chart of the computer program used 
for implimenting equations (1-75) - (1-80) Two additional 
quantities, TP and NUM, are defined. These determine the values 
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of the sampling period T for which the approximations of G w (T) and 

E (T) are calculated. With TP and NUM given the sampling period T 
w 

varies from T = TP to T = TP x Nt)M in NUM increments. 
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1.8 Comouter Program for ComDutation of — Vf ■ and E ^j 

The flow diagram for the computer program that calculates G w 

and E is shown in Fig. 1-4 The program is a straightforward 
w 

implementation of the formulas ( 1-47 ) and (1-63) The subroutine 
STRMAT is. called to calculate the state transition matrix and the 
integral of the state transition matrix over one sampling period. 

Its flow diagram is shown in Fig. 1-5 These matrices are calculated 
by means of the following infinite series. 

i# » AT I 
(j) (T) ** e = 57 + 

■0(T) = J e AX dA = 

J 0 


AT A 2 T 2 

TT+ ST + 


(1-81) 


2 2 3 

IT AT_ A T 

11 21 31 


(1-82) 


Let M be the kth partial sum of the <J> (x) series 


- 

kio k! 

and R the remainder matrix 


(1-83) 


00 _k_k 

R- I 

u kl 


(1-84) 


k=K+l 

Using the upper bounds established by M. L. Liou [1] and W. Everling 
[2] the elements, r\ j , of R may be bounded by 


|r ij 


(AT) 


K 


Kl 


|a| [t 

K+l 


1-e 


(1-85) 


where 


a It 

£ = K+2 


( 1 - 86 ) 
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The error test for terminating the summation of the series (1-81) 
that 



a m. . 
ID 


(1-87) 


where a = 10 6 and m. . are the elements of M. The above inequality 

3-D 

must hold for all i, j except for elements m, . that are close to zero 

13 

If 

| m ij| < 2 - 10- 30 (1-88) 


the test(l-87)£s bypassed. 
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INITIALIZATION 

<f> = 1 
e" = it 

wl = IT 
K = 2 


Test < a W2(i, j) 


N 

N 

ANORM = l 
i=l 

1 1 a ■ • 1 
i=l 13 


W2 = A * Wl 
<J> = <J> + W2 
Wl = W2 * T/K 
~Q = ¥ + Wl 


K > ANORM * T 


N 

N 

BNORM = l 

I | W2 ( i , j) | 

i=l 

j-1 


or 

j = j + 1 


is i = N 
No and 
j = N 


Return 


K > K MAX 


es 


error 


K = K+l 


ANORM * T 


K + 2 


BNORM * ANORM * T 
Test = (K+l) (1-e) 


W2 (i , j) > g 


No 


Figure 1-5 Flow chart for calculation of the state 

transition matrix and its integral, STRMAT 
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1.9 Digital Redesign of the Simplified One-Axis Dynamics of 

the Skvlab Satellite . 

In this section digital redesign with the partial matching metfod is 
applied to the simplified one-axis dynamics of the Skylab Satellite. 

The block diagram of the simplified Skylab Satellite is shown in 
Fig. 1-6 . - The constants of the system are given as 

I = 970,741 Kg - m 2 

3 

K = 11.8 x 10 n - m/rad 

= 125.9 x 10 3 n - m/rad/sec 

The feedback gains g^ and g 2 are to be selected so that the 
damping ratio of the system is 0.707. For zero steady-state error, 
g^ = 1 . The closed- loop transfer function of the system is 


K, 


X (s) 

1 JJ 

iJj (s) 2 

V is +K ig2 s + K 0 


11800 


970741s +125900g 2 s+11800 


0.01215 


s +0.129g s+0. 01215 


(1-89) 


Therefore, the natural undamped frequency is 0.11 rad/sec, and 


g 2 = 1.205, K^g 2 = 151,800. 
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and 


Then 


The state equations of the system are 


X 1 = X 2 


X 2 = 1 ” 


(1-90) 


U " Vc • g lVl “ g 2 K l X 2 


E\p 

r c 

- G 11 (0)x 1 - G 12 

(0)x 2 


(1-91) 

0 

A 

f > 

0 




B = 



(1-92) 

0 

oj 

1/970, 741j 



1 = 

tG U (0) G 12 (0>] 

= [11800 

151,800] 



E (0) = 11800 

(1-93) 

The block diagram for the continuous model of the Skylab Satellite 


is shown in Fig. 1-7 
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The present Skylab Satellite system with the feedback gains 
of Eq.(l-93) has an undamped natural frequency of 0.11 radians/sec. and 
a damping ratio of 0.707. For this system the exact gains G^, E^ 
are computed for T = 1, 2, 3, 4, 5 seconds and H = [0 1] and H = [1 0] 

These gains are listed in Table 1-3. With H = [0 1] the state x^ 

is attempted to be matched and with H = [1 0] the state x^ is 

attempted to be matched. 

The system is simulated on a digital computer for both values of 

H. A sanpling period of T = 2 sec is used for the sampled data system. 

The gains G(0) and E(0) of Eq.(l-93&re used for the continuous system 

and the gains G and E„ (from Table 1-3 for T = 2 sec) are used for 
W W 

the sampled data system. 

Figures 1-8 to l-13show the simulation results for T = 2 sec. and 

H = [0 1] and figures 1-14 to 1-19 show the simulation results for 

T = 2 sec. and H = [ 1 oj . Figures 1-8 and 1-14 show the state 

trajectories x. (t) and x. (t) for the two cases of H respectively, 
lc Is 

Figures 1-9 and 1-15 show the error (x^ c (t) - x^ s (t) ) for the two 
cases of H. The remaining figures show the same results for the 
state x^ and control u. 

The simulation results show that 

i) the states are more closely matched with H = [0 1] than 

with H [1 o] 311(1 

ii) a sampling period of T = 2sec. appears to be adequate for 

digital redesign of the Skylab Satellite system of Fig. 1-6. 
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T (seconds) 

H 

Table 1-3 

W T) 

W T) 

V T) 

1 

[0 

1] 

10901.5 

145840 

10901.5 

2 

M 

10051.2 

139921 

10051.2 

3 

it 

9248.45 

134071 

9248.45 

4 

it 

8492.5 

128315 

8492.5 

5 

it 

7782.34 

122674 

7782.34 

1 

ti 

0] 

11197 

147825 

11197 

o 

it 

10618.1 

143867 

10618.1 

3 

ti 

10063.1 

139937 

10063.1 

4 

it 

9531.78 

136048 

9531.78 

5 

n 

9023.72 

132207 

9023.72 

Also G i;l (0) 

<=12 <°) 
G(0) 

= 11800 
- 151,800 
= 11,800 







Figure 1-7 Block diagram of the skylab satellite with 
state feedback for damping ratio of 0.707. 
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Figure 1-13 Error, u^t) - u g (t), in the control trajectories of the 
simplified one axis dynamics of the Skylab Satellite. 
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2. DIGITAL APPROXIMATION BY POINT-BY-POINT STATE COMPARISON 
WITH EXACT MATCHING AT MULTIPLE SAMPLING PERIODS 

2.1 Introduction 

The problem is that of approximating a continuous-data system of Fig. 2-1 
by a sampled-data model of Fig.)2-2 by comparison of the states of the two 
systems. 

The continuous-data system is described by the following time-invariant 

( 2 - 1 ) 
( 2 - 2 ) 

where a^,(t) denotes an n-vector, ii(t) and ir(t) are m-vectors, E^ and G^ are 

the gains of the forward and the feedback paths, respectively. 

The dynamic equations of the sampled-data model are: 

k (t) = Ax (t) + Bu (t) '(2-3) 

-s ~ s — s 

u (kT) = E(T)r(kT) - G(T)x (kT) - - (2.-4) 

s .s 

for kT _< t _£ (k + 1)T. The sampled-data system is also of the nth order with 

m inputs. The outputs of the sample-and-hold devices are a series of step 

functions which are denoted by the vector u (kT) for kT t < (k + 1)T. 

s 

The objective is to find E(T) and G(T) so that the states of the sampled- 
data model are as close as possible to that of the continuous-data system. 

It was pointed out in Chapter 1 that exact matching of the states of 
the two systems cannot be achieved in general for an n-state and m-input 
system, unless n = m and the following matrix is nonsingular. 


dynamic equations : 

x (t) = Ax (t) + Bu(t) 
— c — c — 

u(t) = E r(t) - G x (t) 



In this chapter we shall show that for a continuous-data system with n 
states and m inputs, the states can be matched exactly, in principle, subject 
to approximations, by the states of a sampled-data system every N sampling 
instants, where 

n/m <_ N •< .n/m + 1 (2-6) 

For instance, for a second-order system with one input, it is possible to match 
the states of the two systems every two sampling instants. On the other hand, 
if the sampling period of the sampled-data system, T, has been fixed, it is 
necessary only to change the forward gain E(T) and the feedback gain G(T) in 
between the sampling instants, at N equal intervals, to obtain exact matching 
of states at the sampling instants. It will also be seen that only the 
values of the states at the sampling instants, kT , k = 1, 2, .... need to be 
fed back to generate the control signals k = 1, 2, ... 
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2.2 Exact Solutions for the Forward and Feedback Gain Matrices 
(Second-order, single-input system) 

In this section we shall derive the optimal gains for the exact matchings 
of states at multiple sampling instants using a second-order system with a 
single input (n = 2, m = 1). 

Let us define 


<KT) - e AT 

(A-BG )T 
<f> c (T) = e C 


(2-7) 

(2-8) 


0 (T) 

c 


(•T (A-BG )A 
e C dAB 

J 0 


(2-9) 


Then, the solution of the state at the sampling instant t = (k + 1)T for the 
continuous-data system of Eqs . (2-1) and (2-2) is 

x [(k + 1)T] = <p (T)x (kT) + 6 (T)E r(kT) (2-10) 

c C C CO 

where it is assumed that r(T) = r(kT) for kT < x < (k + 1)T. This approximation 
is necessary so that r(x) can be factored out of the integral of 0^(1), and 
the solution will not be input dependent. 

Similarly, the state of the continuous-data system at t = (k + 2)T is 
obtained as 


xJU + 2 ) T ] = (jr(T)x c [(k + 1)T] + 0 c (T)E c r[(k + 1)T] 

= <p (2T)x (kT) -4- cp (T) 0 (T)E r(kT) + 9 (T)E r[(k + 1)T] (2-11) 

c — c c c c c c 

For the sarapled-data system, the solution of Eq. (2-3) for the time interval 

kT < t < (k + 1)T, k = 0, 1, 2, ... , is 

x [(k + 1)T] = <p(T)x (kT) + 0(T)u,(kT) (2-12) 

' S w 5a r 
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where u (kT) denotes the output of the new-order hold for the specific time 
s 

interval. For analytical reason we are using the control signal u (kT) in 

s 

the derivation instead of expressions involving the input r(f). 

For two sampling instants, the state vector at t = (k + 2)T is written 

x [ (k + 2)T] = <j>(2T)x (kT) +- <j)(T) 6(T)u (kT) + 0(T)u [ (k + 1)T] (2-13) 
S S 'S s ... 

where u [ (k + 1)T] denotes the constant output of the zero-order hold for the 
s 

time interval kT _< t < (k + 1)T. 

Now let us assume that the states of the second-order continuous-data 
system and that of the sampled-data system can be matched at t = kT and at 
t = (k + 2)T, we let 

x (kT) = x (kT) (2-14) 

s c ■ — “■ 

and 

x [(k + 2)T] - x [(k + 2)T} (2-15) 

- S “C 

Then Eqs. (2-11) and (2-13) give 

<f> (2T)x (kT) + <j> (T)0 (T)E r(kT) + 9 (T)E r[(k + 1)T] 

C — S C C C 0 0 

= <{>(2T)x (kT) + <KT)0(T)u (kT) + 0(T)u [ (k + 1)T] (2-16) 

s s s 

Let 


r(T) = [<J>(T)0(T) 6 (T) ] 

(2-17) 

P(T) = r _1 (T)[<j> c (2T) - cp(2T) ] 

(2-18) 

S. (T) = r -1 (T)4> (T)0 (T)BE 
1 c c c 

(2-19) 

S„ (T) = r _1 (T)9 (T)BE 
l c c 

/ (2-20) 


Notice that T(T) is the "controllability 1 * matrix of the sampled-data system in 
Eq. (2-12) Therefore, if the sampled-data system is controllable with the 
sampling period T, T(T) is nonsingular. 



With the matrices defined in Eq. (2-17) , (2-18) , (2-19), and (2-20) , 

Eq. (2-16) leads to 


u (kT) 
s 


u s [(k + 1)T] 


P(T)x (kT) + S 1 (T)r(kT) 4-S_(T)r[(k + 1)T] 
—s i 2 


( 2 - 21 ) 


where k=2j,j=0, 1, 2, .... 

The last equation is also written as 


u g (kT) = [1 0]{P(T)x s (kT) + S 1 (T)r(kT) + S 2 (T)r[(k + 1)T]} (2-22) 

uj(k + 1)T] = [0 l]{p(T)x s (kT) + S 1 (T)r(kT) + S 2 (T)r[(k + 1)T]} ( 2 


Let 


u (kT) = - G x (kT) + E n r(kT) + F.r[(k + 1)T] (2-24) 

s u— S U u 

where Gg, Eq , and Fg represent the constant gains which are effective during 

the time intervals of 2kT < t < (2k + 1)T, k = 0, 1, 2, ... . Comparing 

Eqs. (2-22) and (2-24) , we have 

G q = -[1 0]P(T) (2425) 

E 0 = [1 0JS x (T> (2-26) 

F 0 - [1 0]S 2 (T) (5-27) 

However, Eq. (2-24) indicates that u (kT) is a function of both r(kT) 

s 

and r[(k + 1)T] , which means that in order to determine the control u (kT) 

' s 

at t = kT, the input r(t) must be known at the next sampling instant 
t ■ (k + 1)T. This implies that either the input signal must be known in 
advance over each interval [kT, (k + 2)T] or a predictor must be constructed 
to predict r[(k + 1)T] on line. 
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In a similar fashion, Eq. (2-23) is written 

u [(k + 1)T] - - G..X [(k + 1)T] + E.r(kT) + F.r[(k + 1)T] (2-28) 

s 1— s - 1 1 

where G^, E^, and represent the constant gains for the time interval 

(2k + 1)T < t < 2(k + 1)T, k « 0, 1, 2 

Substituting Eqs.;(2-24) and (2-12)into the last equation, 

we have, — 

u s [(k + 1)T] = -G 1 [4>(T) - 6(T)G 0 ]x g (kT) 

+ [E 1 - G 0(T)E Q ]r(kT) +' [F ]L - G 1 0(T)F () M(k+l)T] (2-29) 

Comparing Eqs. (2-17) and (2-23) yields 

G 1 = -[0 l]P(T)[(j)(T) - e(T)G 0 3 -1 (2-30); 

E 1 = [0 1]S 1 (T)+ G 1 0(T)E o (2-31) 

F x = [0 1]S 2 (T) + G^0 (T) F q (2-32) 

It can be shown that the condition for the matrix [<f>(T) - 0 (T)Gq] to 
be nonsingular is that the sampled-data system be completely controllable. 

- Note that the control u g [ (k + 1)1] may be expressed in terms of the state x (kT) 
as in Eq. (2-23), or in terms of the state ^[(k+l)!] as in Eq.’(2-28). Thus, the 
state variable can be sampled either every kT or every 2kT, k = 0, 1, 2, ..., fori 
feedback purpose. \ be 

A simpler solution which does not require the use of a predictor can be 

obtained if we assum e th at - 

r [ (k + 1)T] = r (kT) = r(l) ' (2-33) 

for kT £ X < (k + 2)T. Then Eq.^ (2-16) becomes 

<]) c (2T)x s (kT) + 0 c (2T)E c r(kT) 

= <|>(2T)x (kT) + 4>(T)0(T)u (kT) + 0(T)u [ (k + 1)T] (2-34) 

b S S 
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where 

2T (A-BG ) X 

0 (2T) = e C dXB 

C Jo 

Therefore, Eq. (2-34) ± s written 
u g (kT) 

- P(T)x (kT) + S(T)r(kT) 

~ s 

u g [(k + 1)T] 
where k= 2j, j =0, 1, 2, . . . , and 

S (T) = 0 (2T)E 
c c 

Equation (2-36) can be written as 
u g (kT) = - G^CkT) + E Q r(kT) 

u [ (k + 1)T] = - G-X [(k + 1)T] + E. r (kT) 
s i-~s i 

Comparing Eq. (2-36) with Eq. (2-38) yields 
G q - - [1 0]P(T) 

E q = [1 0]S(T) 

Equation (2-39) is rewritten as 

u g [(k + 1)T] = - G^d) - 0(T)G Q ]x g (kT) 

+ [E 1 - G 1 0(T)E o ]r(kT) 

Thus, comparing Eq. (2-36) with Eq. (2-42), we have 
[0 1]P(T) = - G 1 [4>(T) - 0(T)G q ] 

[0 1]S(T) - - G. l 0(T)Eq + E ± 


(2-35) 

[ (2-36) 

(2-37) 

(2-38) 

(2-39) 

(2-40) 
( 2-41 ) 

(2-42) 

( 2 - 43 ) 

(2-44) 


I 
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Solving for G 1 from Eq. (2-43) and E ][ from Eq. (2-44), we get 

G 1 = - (0 l]P(T)[<j>(T) - 9 (T)G q ] _1 (2-45) 

E x - [0 1] S (T) + G^^Eq (2-46) 

Therefore, for the second-order system with single input, if the states 
can be matched every two sampling periods , the gains G(T) and E(T) must 
change twice during each basic period of 2T. Then, during the subsequent 
sampling periods these gains go through the same sequences repeatedly. In 
other words, referring to Eq. (2-4), 

G(T) = G q 2kT < t < (2k + 1)T (2-47) 

G(T) - G ± (2K + 1)T < t < 2(k + 1)T (2-48) 

k-0, 1, 2, ... . In general, if the states of the continuous-data system 
are exactly matched with those of the sampled— data system every NT seconds, 
where N is an integer greater than unity, the feedback gain would have to 
assume N different values, G^ , G^, ..., G^_^, over the basic period, NT, 
with changes occuring at the sampling instants T. Furthermore , 


G k =G k+N < 2 - 49 > 

for k = 0 , 1, 2, ... . Similar properties may be stated for the forward gain 

E(T) of Eq. (2-4) 

As an alternative to sampling and feeding back the states at. every 
sampling period T, we can assume that the gain matrices are changed every T 
but the states are sampled only at t = 2kT. Then, in addition to Eq. (2-38) 
Eq. (2-39) is written as 


u g [(k + 1)T] = - G 1 x g (kT) + E x r(kT) 


(2-50) 
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where k=2j,j = 0,l,2, ... . The gains G Q and E Q in Eq. (2-38) are 
still given by Eqs. (2-40) and (2-41) , respectively, but comparing Eq. (2-50) 
with Eq. (2-36), we have 

G - - [0 1]P (T) (2-51) 

and 

E x - [0 1]S(T) (2-52) 

It should be noted that feeding back the states every T seconds and 

feeding back the states every 2T would yield the same responses for the system 

at the sampling instants, although the control signals and the gains for the 

two schemes are different. However, using x (kT) k = 0, 2, 4, ... means that 

— s 

the states may be sampled only half as frequently in the second-order system, 
but the gain matrices must be changed during the sampling period. In practice 
the choice of either method may depend upon the applicability of the methods. 

For instance, if the sampling period T is fixed a priori, then we may say that 
for the second-order system considered it is possible to match the states at 
every other sampling instants only, by feeding back the state variables and 
changing the gains at each sampling instant. However, in order to match the 
states at every sampling instant we can change the forward and feedback gain 
matrices once in between the sampling instants. If we use Eqs. (2-51) and 
(2-52) for this purpose no additional samplings of the states are necessary. 

Although it may be more convenient to sample the states every N sampling 
periods and change the feedback and forward gains every sampling period to 
match the states every N sampling periods, it should be noted that the system 
is operating open loop over these N sampling periods. If the states and gains 
are changed at the end of every sampling period the system will only be operating 
open loop over a single sampling period. This second method requires that the 
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states be sampled more frequently but, in general, it should be less 
susceptible to noise perturbations and variations of system parameters. 

Another interesting observation is that the digital point-by-point 
state matching problem considered here is closely related to the time- 
optimal control of sampled-data control systems. For simplicity of illustration, 
we may consider the regulator problem by setting r(kT) = 0 for all k. Then, 
the optimal control for point-by-point state matching for the second-order 
system is obtained from Eq. (2-21) 


K<kT) 


= P(T)x (kT) 


u s [(k + 1)T] 


» [<j>(T) 9 (T) 9 (T) ] _1 [<fr (2T) - <|>(2T)]x (kT) (2-53) 

c —s 

In sampled-data systems theory, it is well known that given a second- 
order system with initial state x (kT) , the optimal control to bring x (kT) 
to the equilibrium state () in minimum time is 


= [<KT)0(T) 0(T)]" 1 x s (kT) (2- 54 ) 

Therefore, Eq. (2-47) can be interpreted as the optimal control which brings 

the difference of <f> (2T)x (kT) and 4>(2T)x (kT) to zero in minimum time. For a 

c — c — s 

second-order system, if the amplitudes of the controls at various sampling 
instants are unconstrained, the mimimum (and maximum) time to achieve state 
matching is two. However, if the. controls are subject to amplitude constraint, 
more sampling periods are required for matching. 
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2 . 3 Exact Solutions for the Forward and Feedback Gain Matrices 
(General Case ) 

Now consider that the continuous-data system of Eqs.(2-1) and (2-2) 
has n states and m inputs. Then, if the sampled-data system of Eqs . (2-3) 
and (2-4) is controllable for the sampling period T, the states of the two 
systems may be exactly matched every N sampling instants; where N is given by 


Eq. (2-6) . 

The states of the continuous-data system evaluated at t = (k + N)T are 
expressed as 


x t(k + N)T] = <f> (NT)x (kT) + 6 (NT)E r(kT) 
— c C — c c c— 


(2-55) 


where 


<f> c (NT) = e 


N(A-BG )T 
c 


e c (NT) = 


NT (A-BG )A 
e C dXB 

0 


(2-56) 


and 


By recursion, the solution of Eq. (2-16) at (k + N)T is 

N-l : 

x [(k + N)T] = <j> (NT)x (kT) 4- £ <p [ (N - i - l)T]6(T)u [ (k + i)T] (2-57) 


— s 


i=0 


Now equating Eq. (2-55) to Eq. (2-57) and letting x^kT) = Xg(kT), we have 

N-i 

[<{> c (NT) - (|)(NT)3x s (kT) + 6 c (NT)E c r(kT) = £ <f> [ (N - i - l)T]0(T)u s [(k + i)T] (2 _ 58) 


The last equation can be written as 
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[<f>[(N - 1)T]0(T) «j)[(N - 2)T] 6(T) . . . <fr(T)0(T) 6(T>] 

n x Nm ► 

1 

< ■■■ Nm x 1 *■ 

= [<J> c (NT) - <{)(NT) Jx^ (kT) + 0 c (NT)E ( jr(kT) , (2-59) 

Let - .... 

T(T) = [<(> [ (N - 1)T]0(T) <J>[(N - 2)T]0(T) . . . (j)(T)0(T) 0 (T) ] (2-60) 

which is not a square matrix unless n = Nm. However, if the sampled-data 
system is completely controllable, F(T) must be of rank n. Let A(T) be an 
n x n matrix which is formed by use of n linearly independent columns of T(T). 

In general, there may be a nonuniqueness in the selection of A(T) if NM > n and 
there are more than n independent columns in T(T). Now let IJ represent at 
n x 1 vector which contains the elements of the control vector of Eq. (2-59) 
that correspond to the columns of A(T) taken from T(T), which means that 
Nm - n elements of the control vector can be set arbitrarily. Then, Eq. (2-59) 
leads to 

_U= A _1 (T) [<J> c (NT) - (KNTJlXgCkT) + A _1 (T)0 c (NT)E c r (kT) (2-61) 

The procedure of expressing the elements of _U in terms of state feedback in 
the form of the following equations is conceptually simple: 

u [ (k + j)T] = - G x (kT) + Er(kT) (2-62) 

^~S J” S J 


V^CkT) 

u [(k + 1)T] 

5 

u [ (k + N - 2)T] 
— s 

u [(k + N - 1)T] 
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or 


u [ (k + j)T] = - G.x [ (k + j)Tj + Er(kT) 

—3 j-S j— 


(2-63) 


j = 0, 1, 2, .... N - 1. 

However, because of the nonuniqueness in the solution general closed form 
solutions for the gains are dif ficult^to express. 

Let us consider the case when n/m is an integer; then N = n/m. In this 
case, r(T) is an n x n square matrix, and the N optimal control vectors, 
u^dcT), u^[(k + 1)T], •••» Ug[(k + N - 1)T] , are uniquely determined from 

h (kT) 

— s 

u [(k + 1)T] 


n x l 


u [(k + N - 2)T] 

5 

u [ (k + N - 1)T] 
— s 


= P(T)x (kT) + S(T)r(kT) 


(2-64) 


where 


P(T) = r _1 (T)[(|) c (NT) - <f>(NT)] 

S (T) = r _1 (T)0 (NT) E 
c c 


(2-65) 

( 2 - 66 ) 


Then the control vector, u^Kk + j)T), j = 0, 1, 2, ...» N - 1, can be 

A A 

expressed by Eq. (2-62)or Eq. (2-63) The matrices , G^ , G^. , , and , 

represent the gain matrices which must change N times during each sampling period 
NT. During subsequent matching periods of NT to 2NT, 2NT to 3NT, ..., etc., 
these gain matrices go through the same N sequences of variations for each 
period. 


Let us define the m x n matrix 1^ as 
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l 


• 0 . I • 1 . I • 0 . 

• • • ■ • • • 

0 0 | o 1 I 0 0 

■*- m(j + l)->- m -*• — 2m- 

columns columns - columns 


m 


(2-67) 


3 »!»••• 


then the control vectors of Eq. (2-64) can be written as 

u [(k + j)T] = I,[P(T)x (kT) + S(T)r(kT) ] 

— s . j — s — 

Now comparing Eq. (2-62) with Eq. (2-68), W e have 


( 2 - 68 ) 


and 


G = - I.P(T) 


E^ = I S(T) 


(2-69) 


( 2 - 70 ) 


for j = 0, 1, 2, . . . , N - 1, 

which are relatively simple solutions for the gain matrices. 

However, if the control vector is defined by Eq. (2-63), the solution is more 
complex. 

For j = 0, comparison of Eqs. (2-63) and (2-68) gives 

« 0 ■ - I 0 ?CT) (2-71) 

E 0 " I 0 S(T) (2-72) 

As expected, when the exact matching of states can be accomplished in 

one sampling period, the solutions from Eqs. (2-62) and (2-63) are identical. 

For j = 1, Eq. (2-63)gives 

u [(k + 1)T] = - G. [<j>(T) - 0(T)G n ]x (kT) 

— S 1 u — s 


+[E 1 - G 1 E Q ]r(kT) 


(2-73) 


I 
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Comparing Eq. (2-73) with Eq. (2-68), with j - 1, we have 


G 1 = -1^(1) [cKT) - e(T)G 0 3 -1 (2-74) 


E 1 = I 1 S(T) + (2-75) 

For the general case. 


j 

x s [(k+ j)T] = n [<KT) - 8(T)G j _ i ]x s (kT) 
j" 1 j-P 

----- - + l n [4>(T> - 0(T)G ] 0(T)E r(kT) 

p=l i=l J P 


+ 0(T)E_j_ 1 £(kT) .(2-76) 


Substituting Eq. (2-76) into Eq.(2-63), we have 

j 

ju^CCk + j)T] = - Gj 11^ [4>(T) - 0(T)G J _ i ]x s (kT) 

j-l j-P 

- G l n [<KT) - 0(T)G. ]0(T)E r(kT) 

J p=l i=l 3 1 P 


- G.0(T)E ,r(kT) + E.r(kT) (2-77) 

j j-1— J- 

Comparing Eq. (2-77) with Eq. (2-68) the gain matrices are obtained as 


G = - I ,P (T) { n [<P(T) - 0 (T)G ,] } -1 (2-78) 

j j i= l 3 ~ i 

j-P 

E, = I S(T) + G, l n [4>(T) - 0 (T)G ] 0 (T)E 

J J J p-1 i=l ^ _i P_1 


+ GjQ (2-79) 


j = 0, 1,2, ••>, 


N - 1. 
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2.4 The General Gamma Matrix 

In section 2.3 we considered the case when n = Nm and the system was 
controllable in N sampling periods. A more general case is when n/m is 
not an integer or when n = Nm but the system is not controllable in N 
sampling periods. In both these cases N must be chosen such that Nm > n. 

The matrix T(T) given by Eq. (2-60) is then nonsquare. As mentioned previously 
if an integer N can be found such that the system is controllable in N 
sampling periods and Nm > n then a nonsingular matrix A(T) of dimension n x n 
can be formed by selecting n linearly independent columns of T(T). Equation 
(2-59) is then written as 

A(T)U(T) + A(T)U(T) = [<() (NT) - 4> (NT ) ]x (kT) + 0 (NT)E r(kT) ( 2 -80) 

c s c c 

where u(T) is an n-dimensional control vector whose elements correspond to 
the n linearly independent columns of T(T) that were selected to form A(T) . 

The remaining Nm - n columns of T(T) form the n x (Nm - n) matrix A(T) , 

and the remaining Nm - n controls form the corresponding control vector U(T). 

Since the controls _U(T) may be set arbitrarily, let U(T) = 0. Solving Eq. (2-80) 

gives I _ . 

U(T) = P (T)Xg (kT) + S(T)r(kT) (2-81) 

where 

P (T) = A' 1 (T)[<|) c (NT) - 4) (NT)] 

S(T) - A _1 (T)0 c (NT)E c 

First consider the case when the control is expressed as 

A A 

^[(k + j)T] = - GjX g (kT) + E^r (kT) (2-82) 


j ° 0, 1,2, ..., 


N - 1. 
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It is sufficient to consider u (t) over the first sampling period since 

s 

A A 

the gains G^ and E^ repeat themselves every sampling period. 

In section III we equated the first m columns of _U(T) with u (0) to 

"" s 

A A 

find Gq and E^. However here the vector JLJ (T) in Eq. (4-2) may not contain 

all the controls for the first sampling period, u, (0) i = 1, 2, ... m. 

^ — i 

Some of these controls are contained in the U(T) vector and some the U(T) vector. 


The ones that are in the U(T) vector have been set to zero. In order to 
compare the controls in equations (4-2) and (4-3) it is necessary to add 
back to the U(T) vector the controls that have been set to zero. That is, we 
form the Nm - vector 


U(T) - 


fug(kT) 

u [(k + 1)T] 
s 

u [ (k + N - 1)T] 


(2-83) 


by adding back the Nm - n controls that formed the U vector. We also construct 
an Nm x n matrix P(T) by inserting rows of zero corresponding to the 
additional Nm - n controls. Similarly we construct an Nm x m matrix S(T). Then 


U(T) = P(T)x (kT) + S(T)r(kT) 

— — s — 

The control vectors of Eq.(2-83) can then be written as 

u [(k + j)T] = I [P(T)x (kT) + S(T)r(kT)] (2-84) 

Comparing Eq. (4-3) and Eq. (2-84) we have 

ej-V (T) (2-85 

E j = I j^( T ) (2-86) 

for j = 0, 1, 2, .... (N - 1). 
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Next we consider the case when the control is expressed as 


u [ (k + j)T] = - G.x [(k + j)T] + Er(kT) 
s j s J 


(2-87) 


for j - 0, 1, 2, , (n - 1). 

Equating the control given by Eq. (2-77) and the control given by Eq. (2-84) 
gives expressions identical to Eq. (2-78) and Eq. (2-79) except that the matrices 
P(T) and S(T) are replaced by P(T) and S(T). 

In order to illustrate the matrix manipulations discussed in this section 
we will consider the case of a third-order system having two inputs. 

Equation Eq. (2-59) is written as 


X 11 X 12 X 13 X 14 


X 21 *22 X 23 X 24 

X 31 X 31 X 33 X 34 



^(0)' 


u 2 (0) 


u x (T) 


u 2 (T), 


- [<f> (NT) - 4> (NT) ]x (kT) +0 (NT)E r(kT) (2-88) 
c s c c 


Since n = 3 and m = 2 the T(T) matrix is nonsquare. If the first three columns 
of T(T) constitute a nonsingular matrix, we can define A(T) as 


Also 




X 12 

X 13 

A(T) = 

X 21 

X 22 

X 23 


X 31 

V 

X 32 

X 33 


U’(T) = [UjW) u ? (0) -u^CT)] 

Setting u 2 (T) = 0, Eq. (2-88) is written as 


..(2-89). 


/ (2-90) 


A(T)U(T) = [ cj) c (NT) - (}) (NT) (kT) + 0 c (NT)E c r (kT) 


(2-91) 



where 


P(T) - A -1 (T) [<J> (NT) - 4> (NT)] 
c 

S (T) = A _1 (T)0 (NT)E 
c c 

By adding a row of zeros to the P(T) and S(T) matrices Eq. (2-92) may 
written for the complete control vector for both sampling periods. 


u x (kT) 

u 2 (kT) 

u x [(k + 1)T] 
u 2 [(k + 1)T] 


= P (T)x (kT) + S(T)r(kT) 
— s — 


where 


P(T) = 


P 11 

P 12 

P 13 


b ll 

S 12 

P 21 

P 22 

p 

23 


S 21 

S 22 

P 31 

p 

32 

P 33 

S(T) = 

S 31 

S 32 

_ 0 

0 

o . 


. o 

0 

* 


The G , E, and G. , E. gains may then be found by using Eqs. (2-69), (2 

J 3 J • J — - ---■ 

and (2-72)> (2-73) respectively. 


(2-93) 

(2-94) 

be 


-70) 
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2 . 5 Simulation of Systems 

Digital Redesign of the Simplified One-Axis Dynamics of the Skylab Satellite 
The dynamic equations of the simplified one-axis model of the skylab 
satellite were discussed in detail in Chapter 1. 

The state equations for the continuous system are 

dx 

-7= - [A - BG ] x ( t ) + BE r (t) 
dt c — c 

where 


r(t) = unit step function. 



r o 

1 


’ 0 ' 

A = 



B = 



0 

1 


1/970,741 


G - [11800 151,800] 

c 

E = 11800 
c 

The method of multiple sampling periods discussed in sections 2.2 - 2.4 was used 
for the digital redesign of the above system. Since the system is of the 
second order with one control, the responses of the continuous and digital systems 
can at best be matched at the end of two sampling periods. For the skylab 
example the digital system was controllable in two sampling periods so that exact 
matching could take place at the end of the second period. The sampling 
period was chosen as T = 1 sec and N = 2. Therefore, matching occurs at 
t = 2, 4, 6, ... sec. The feedback gains and the forward gains are 

. G Q - [11185. 147812.] 

G = [10639.6 144149.] 

E q = [11185.] 


E - [10639.6] 



Gq and E q are used over the sampling periods kT for k = 0, 2, 4, ...» and 
and are used over the sampling periods kT for k = 1, 3, 5 • •• • 

The simulation results are shown in Figs. 2-4 ' through ; 2 _g ! f° r the states 
and controls and the errors between the states of the digital and continuous 
systems. The matching of the states can best be seen in Fig. 2-7 At the 
end of the first sampling period the error is negative, however, at the end of 
the second sampling period the error is close to zero. 

Table 3-1 in the next chapter compares the method of partial matching to the 
method of multiple sampling periods by giving the maximum error between the states of 
the digital and continuous systems at the end of one sampling period for the 
method of partial matching and at the end of two sampling periods for the 
method of multiple sampling periods. 
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3. DIGITAL APPROXIMATION BY POINT-BY-POINT STATE 
COMPARISON WITH HIGHER ORDER HOLDS 

3 . 1 Introduction 

The problem is that of approximating a continuous-data system of Fig. 3-1 
by a sampled-data model of Fig. 3-2 by comparison of the states of the systems. 

The continuous-data system is described by the following time-invariant 
dynamic equations: 

x (t) = Ax (t) + 3u(t) (3-1) 

-~c — c ~ 

u(t) = E^r (t) - G^Ct) (3-2) 

where x^ (t) denotes an n-vector, _u(t) and r/t) are m-vectors , E and G are the 
c c c 

gain matrices of the forward and the feedback paths, respectively. 

The dynamic equations of the sampled-data system which is to replace the 
continuous-data model are 

x (t) = Ax (t) + Bu (t) (3-3) 

c S — S ... — “ 

where u^t) is an m-vector and is defined as the output of an (N-l) st-order 
hold, N-l, 2, ... . 

It will be shown that it is possible to match all the states of the system 
in one sampling period if a hold device of sufficiently high order is utilized. 

The order of the hold which is necessary to match the states, depends upon the 
ratio of the number of states and the number of controls of the system. For 
the system considered which is of the nth order with m inputs, it is possible 
to match all the states at 'every sampling instant with at least an (N - l)st-order 
hold, where n/m <_ N < n/m + 1. 
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Let us consider, for the time being, that an (N - l)st-order hold is 
adequate for exact matching of states at every sampling instant. (The case 
when a higher than (N - l)st-order hold is necessary will be discussed later.) 
The control input of the sampled-data system with an (N - l)st-order hold is 
expressed as 


UgCt) = [E Q + (t - kT)E 1 + 


N-l 


(t - kT) 

(N - 1)! 


E N _l ] ^ (kT) 


- [G q + (t - kT)G 1 + 


(t - kT) 


N-l 


(N-l)! 


S-l^ 


(kT) 


(.3-4) 


for kT _< t < (k + 1)T, where , E , ...» E^ ^ are m x m coefficient matrices, 

and G„ , G, , ..., G„ , are m x n coefficient matrices. 

0 1 N-l 

Let us define the following matrices: 




(mN x m) 


(3-5) 



(mN x n) 


(3-6) 


G 


N-lj 


Then, Eq. (3-4) is written 


u (t) = [I 
— s m 


(t - kT) I 


m 


(t - kT) N_1 

( (N"?V I J [ ^ (kT) 


(3-7) 


- Gx (kT)] 


kT < t < (k + 1)T 
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where I is the indentity matrix of dimension ra x m. 
Now let us define the m x mN matrix 


L (t - kT) = [I (t - kT) I 
mN m m 


(t - kT) N_1 , 
(N - 1)! m J 


Equation (3-7) is written as 


u 


(t) = E(t - kT)r(kT) - G(t - kT)x (kT) 


where 


E(t - kT) = L A1 (t - kT)E 
mN — 

G(t - kT) = L^Ct - kT)G 


(3-8) 


(3-9) 


(3-10) 

(3-11) 


lcT < t < (k + 1)T. 
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3.2 Derivation of the Optimal Gains 

The necessary conditions for the exact matching of all the states of the 

continuous-data and the discrete-data systems at t = (k + 1)T, k = 0, 1, 2, ..., 

for arbitrary inputs r(t) and initial conditions, x (kT) = x (kT) , are 

— — c ~S 


given by Eqs. (1-11) and (1-12), 


(A-BG )T 
c 
e 


AT 

e 


(k+l)T 

JkT 


e A(kT+T-x) 


BG(t - kT) dT 


(3-12) 


where G^ denotes the feedback gain matrix of the continuous-data system, and 
G(t - kT) is given by Eq. (3-11) 


(k+l)T (A-BG ) (kT+T-T ) f(k+l)T 

e C BE dr = 

kT C * kT 


e A(kT+T t ) be(t _ kX ) dx (3-13) 


where denotes the forward gain matrix of the continuous-data system, and 
E(t - kT) is given by Eq. (3-10), Also, in arriving at Eq.(3-13) it has been 
assumed that _r(t) = i:(kT) for kT _< t < (k+l)T. 

If we let A= (k + 1)T - T, Eqs. (3-12) and ( 3-13 ) become 


where 


AT 

e - e 

6 (T) = 
c 

A 

A = A - 

8 (T) = 
c 


AT = e^BGtf - A)dA 


e AA BE(T - A)dA 


BG 


f T 


e^BE dA 

C' 


(3-14) 

(3-15) 

(3-16) 

(3-17) 


and 
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Substituting Eq. (3-11) into Eq. (3-14) we have 


AT AT 
e - e 


fT 


e AX BL (T - A)GdA 
mN — 


(3-18) 


V- ..j 


Similarly, substituting Eq. (3-10) into Eq. (3-15)gives 


9 (T) = | e AX BL (T - A) Ed A 

>0 mN 


(3-19) 


Equation (3-18) represents n equations in nmN unknown variables in (J, 

2 

and Eq. (3-19) represents nm equations in m N unknown variables in E_. From 

2 

the assumption that n/m _< N < n/m + 1, there are at least n variables in 
Eq. (3-18) and at least nm variables in Eq.(3-19). Thus, it is possible under 
certain conditions to solve for G and _E from these equations. 

If we define 


(T) . f T dX 

J 0 . ••• 

i=0 , . .., N - 1, then Eqs. (3-18) and '-( 3-19 )become 


(3-20) 


e AT - e AT = [q 0 (T) q^T) ... q^COlG 
= Q(T)G 


9 (T) = Q(T)E 
c ~ 


,(3-21) 

(3-22) 


Note that 


q (T) = I e AA BdA = 0(T) 
J 0 


(3-23) 


The matrix q^(T), i = 0, 1, ..., N-l , is of dimension n x m, and Q(T) is of 


dimension n x Nm. 
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The expression for q^(T) in Eq. (3-20) can be simplified when i >_ 1 , 


v = (t - X) /i! 


, AO ^ 
dw - e do 


then 


, (T - X) 1 1 „ 

dv (i - ijT dx 


A0 , 

w = I e da 


where, in general, a is an arbitrary constant. 

Substituting the last four equations into eq. (3-20)> and applying 
integration by parts, we get 


f X Act 
e da 

A=T 
B + 

r 

ff "J 

e da 

a 

A=0 J 

0 

L 
i a 


(T - A) 


i-1 


(i - 1) 


-dAB 


Or, 


- -TT 


O 

[ f° AO 

£ 

— -/ 

e daB + 

e do + 

! e da 

J a J 0 

l j a 

J o J 


(T - A) 


i-1 


(i-1)! 


-dAB 


T 1 f° A (J^ , T 1 f° AO , f T [ X A0, (T - A) 1 " 1 „„ 
e dCB + - 77 * e dOB + e da — — r r dAB 

J a J 0 J 0 


l! 


li 


T ^ AO (T- A)'*'"'^' 

e^d dAB 

o JQ 


Now using Eq . ( 3-23 ) , the last equation is written 


q^D = 


0(A) 


(T - A) 
(i-1)! 


i-1 


dA 


. Let 
(3-24) 
(3-25) 

(3-26) 

(3-27) 

(3-28) 


(3-29) 


(3-30) 
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Thus, Q(T) is given by 


Q(T) = [q 0 (T)j q]L (T) : . . . q^CT)] 


= [0(t) ! 


AA 


fT 


rT 


e (x)dx: 


o 


6(A) (T - A)dA I 


* JO 


T (T - A) N “ 2 

e(A) U -AT dA] 


(N - 2) 


(3-31) 


If 0(A) and e are available in explicit form, Eqs . (3-30) and (3-31) are 
useful for evaluating Q(T). However, if 0(A) and e^ have to be calculated 
numerically using an infinite series, an alternate form for q^(x) and Q(T) is 
more useful. This alternate form is now derived as follows. 

The infinite series expression for e^ is written 


AA r aW 

e = l 

j-0 3 

Substituting Eq. (3-32) into Eq. .( 3-20 )yields 


(3-32) 


qi (T) = 


w 

l (T - A) 1 


j=0 j! 


i! 


BdA 


" I 

j=0 


j rT A U (T - A) . 




BdA 


0 


(3-33) 


In order to carry out the integration in the last equation, for i ^ 1, we 
choose 


then 


V « (T - A) /i! 


dw = -rr dA 


dv = - 


(T - A) 


i-1 


(i - 1) 


dA 


J +1 


w = 


(3+1) ! 


(3-34) 

(3-35) 

(3-36) 

(3-37) 
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Substituting the last four equations into Eq . (3-33), and carrying out 
the integration by parts , we have 


m _ y A J A j+1 (T - A) 1 

- l a (J * D! il 


a=t 

B + 

A=0 3=' 


“ . (T ,j+l / T , \1“1 


. I 


i=0 


T , j+1 


, i-1 


0 


(T - A)' 

(j + 1)! (i-1)! dAB 


( 3-38) 


For i 2 , the above process can be repeated with 
.i-1 


v = 


and 


to give 


dw = A 


(T - A)- 

(i - D! 
3+1 


(j+1)! 


, y A j f T Xj+2 - A) 1 

q i <T) ■ /. 0 Jo <J + 2 > ! tt- 


i-2 


2)1 


dXB 


(3-39) 


(3-40) 


(3-41) 


Generalizing the above described process, it can be seen that the 
expression for q^(T) in Eq. (3-33) needs to be integrated by parts a total 
of i times, to yield the result 


00 • r T • i • 

*3 I ,3 +1 


q ± <« = I d^B 

j=0 0 (i + j)! *** 


= j aV +1+1 


i=0 


(i + j + 1)! 


(3-42) 


Thus , 
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00 j j+1 

fjVur B ' 3 -- 3 > 

V T) ■ jo TTTTjT B ^ 


* I 

lN J- 4 „n / • 


j T j+N 


and 


j=0 (j + N)! 

f 


Q(T) = l A j T j+1 I t TtJtT 1 : - 

- |.(j+l)! n * U +2 ) ! n - 


j=0 


(3-45) 


(j+N) ! n 


(3 


where I is the n x n identity matrix, 
n 

The above expressions for Q(T) can be readily implemented on a digital 
computer. 

The desired solutions for G_ and IS are to be obtained from Eqs . (3-21) 
and (3-22) , respectively. To solve these equations, two cases have to be 
considered: 

Case 1. n/m = N. 

In this case Q(T) is of dimension n x n and (3 is also of dimension n x n 
If [Q (T) ] ^ exists, then the desired solutions of Eqs. (3-21) and (3-22) are 

A 

G “ [Q(T)] _1 [e AT - e AT ] (3-47) 

and 

E= [Q(T)] -1 0 c (T) (3-48) 

where E and 6 (T) are n x m matrices. Substituting these solutions into 
c 

Eqs. (3-10) and (3-11) yields the optimal feedback and forward gains 


-46) 
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G(T) = L mN (T)[Q(T)] " 1CeAT '" eAT] 

E(T) = L mN (x)[Q(T)] " 1 c (T) (3 " 49) 

where X = t - kT , kT < t < (k + 1)T , 0 < X < T. 

Case 2 . n/m < N < n/m + 1. 

In this case Q(T) is of dimension n x Km, G and _E are of dimensions Nm x n 
and Nm x m, respectively. Since Nm > n and Nm > m, Eqs . (3-21) and (3-22) have 
more variables then the number of scalar equations. 

Let 

Nm - n = k (3-50) 

Then, the number of excess variables in Eqs. ( 3— 21 )i and (3-22) are kn and km, 
respectively, and the number of columns to be deleted from Q(T) to yield an 
n x n matrix is k. If the rank of Q(T) is n, then k dependent columns of Q(T) 
may be deleted and the corresponding k rows of and 15 may be chosen arbitrarily. 
Let the n x Nm-dimensional Q(T) be written as 

Q(T) = [Q n (T> : Q k (T)] (3-51) 

where 

Q^(T) = n x n nonsingular matrix 
Q^(T) = n x k matrix 

and k is given by Eq. (3-50) In general, given Q(T) with rank n, there is a 
degree of arbitrariness in the selection of Q n (T) which is nonsingular and is 
n x n. 

Let G be the n x n matrix which contains the n rows of G that correspond 
— n 

to the columns of Q (T) , Also G, is the k x n matrix which contains the remaining 

n — k 

rows of G that correspond to the columns of Q^(T). Then, Eq. (3-21) is written 



(3-52) 
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e AT - e AT = Q(T)G 


= tQ n (T) j Q k (T)] 


r G ' 
— n 


£ j 


Similarly, let E and E, be the matrices which contain the rows of E that 
— n — k — 

correspond to the columns of Q (T) and Q^(T), respectively; E^ is n x n and E^ 
is k x m. Then, Eq. (3-22) is written 


e c (T) = Q(T)E 


- [Q (T) : Q, (T) ] 
n • k 


'e ^ 


i4j 


(3-53) 


Since the nk elements of may be chosen arbitrarily, Eq. (3-52) is written as 


e AT _ e AT _ q ( T )n = Q (T)G (3-54) 

k n — n 

Similarly, the mk elements of _E^ are chosen arbitrarily; thus Eq. (3-53) becomes 

e c (T) - Q^T)^ = Q n (T)4 (3-55) 

Since it is assumed that Q (T) is nonsingular, G and E are solved from the 

n — n — n 


last two equations, 


G = [Q (T)] 1 [e AT - e AT - Q (T)G ] 
— n n k — k 


(3-56) 


^ - lQ n (T)]' 1 [e c (T) - Q k (T)E fc ] (3-57) 

Now using Eqs. (3-10) and (3-11) the optimal feedback and forward gain matrices 


are written 


G < T > - - W T) 


r G ' 
~n 


c k 


- L mN (T) 


f [Q n (T)] _1 [e AT - e AT - Q k (T)G k ]] 


^k 


(3-58) 
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E<T) - " W T > 


E 

— n 

4 

\T*-J 


L mN (T) 


tqmrhejT) 

n c 

4 


V T> 4 ] I 


(3-59) 


where T = t - kT, kT £ t < (k + 1)T, 0 < T < T. 

In general, if there are no amplitude constraints on the control inputs, 
we may set the elements of and arbitrarily to zero. Then, Eqs. (3-58) 
and (3-59) become 


f 


G(T) ’ L „N (t) 


[Q n (T)T 1 [e AT - e AT ] 


0 


E < T > - W T) 


[Q (T)] -1 0 (T) 
n c 


(3-60) 


(3-61) 


0 
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3 . 3 Example of Digital Redesign with the Method of Point by Point State 
Comparison by use of Higher-Order Holds . 

In this chapter the simplified one-axis dynamics of the Skylab Satellite 
System is digitally redesigned by the point -by-point method using higher 
than zero-order holds. The optimal gain matrices G(t) and E(t) are 
Computed, and the system responses are simulated with unit-step inputs. 


This system has been simulated earlier and its system and 
control matrices are 



'0 1 


' ■ 

0 

A = 

o 0 

, B = 

1 



[970,741 


(3-62) 


G c = [11800 151,800], E c = [11800] 


The continuous system gains are 

(3-63) 

Since n « 2 and m = 1 a first-order hold should be adequate to match 
the states at every sampling instant. The feedback and forward gains of 
the digital system are defined by use of Eqs. (3-4) and (3-9) , 

(3-64) 


G(t - kT) = C4 q + (t - kT)G 
E(t - kT) = E q + (t - kT)E 1 
kT < t < (k + 1)T 


Thus 


V 


V 

? E. - 

G, 


E, 

l 1, 


l lJ 


G = 


The state transition matrix is 


(3-65) 


(3-66) 


AT 


1 T 
0 1 


(3-67) 
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and from Eq. (2-23) 


q 0 (T) - 0(T) = 


(T 


0 


e A *Bd* = 1 


Also, from Eq. (2-30) 


q x (T) = 


rT 

0 


9 70,741 


T 

6 

T 2 


f T 2 

2 

IX 


q 0 (X)dA 970,741 
Using Eq. (2-31), we have 


Q(T) [q Q (T) q x (T)] 970,741 


T _ 
2 


Thus , 


[Q(T )]' 1 = 970,741 


6_ z 2 

t 2 T 

12 6 _ 

^3 „2 

T . T 


r 

6 

2 


(3-68) 


(3-69) 


(3-70) 


(3-71) 


The solutions for (2 and _E are obtained from Eqs . ( 2 - 47 ) and (2-48)» 
respectively ; 


G = [Q(T)] _ 1 [e AT - e AT ] 


T 

-12 


-_2 

T 


r AT 
[e 


AT. 
e ] 


(3-72) 


= 970,741 
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E = [Q(T)] _ 1 e c (T) 



= 970,741 




For T = 2, we have 


11752 

151758' 


M 

-1700.7 

-11837 


i G i. 


(3-73) 


( 3-74 ) 


11752 


f E °l 

r 1700.7 




Thus , 

G Q = [11752 151758] 

G x = [-1700.7 -11837] 

E q = [11752] 

E 1 = [-1700.7] 


(3-75) 


1(3-76) 


The optimal gains are obtained from Eqs . (3-64) and '(3-68) 

G(t - kT) = [11752 - 1700.7 (t - kT) 151758 - 11837 (t - kT) ] 

:( 3-77 ) 

E(t - kT) = [11752 - 1700.7 (t - kT) ] 
kT < t < (k + 1)T 
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The continuous and digital data systems are simulated on a digital 
computer. A unit-step input is applied in both cases and a sampling period 
of T = 2 sec. is used for the sampled-data system. The results of the 
simulation are shown in Figs. 3-3 through 3-8 These include the 
state trajectories, control trajectories and the errors in the state and control 
traj ectories. 

Comparing these simulation results with those of the partial- 
matching method of Chapter 1, it can be seen that the 
error between the states and the controls of the continuous and sampled- 
data systems is lower with this method. Theoretically, the 
errors should be zero at each sampling period but due to the numerical 
processes of simulation a small amount of error is still present. Table 3-1 
shows the maximum errors at the sampling instants and in between the sampling 
instants for the exact as well as the partial matching cases. 



00 *s 
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Unit -Step Input 
T = 2 sec. 

G(t - kT) = [11751.9 - 1700.7 (t - kT) 151758 - 11837 (t - kT) ] 

E(t - kT) = [11751.9 - 1700.7 (t - kT) ] 



Time (Sec.) 

Figure 3-7 . Simplified Skylab Satellite 

Comparison of Controls , u and u • 




Unit-Step Input 
T = 2 sec. 

G(t - kT) = [11751.9 - 1700.7 (t - kT) 151758 - 11837 (t - kT) ] 

E(t - kT) = [11751.9 - 1700.7 (t - kT) ] 



Figure 3-8 


Time (Sec.) 


Simplified Skylab Satellite. 
Error, u - u . 



Table 3.1. Comparison of errors between the states of the continuous and 
digital systems for the Simplified Skylab Satellite 


Method of 
Digital Redesign 

Maximum Error at the Sampling Instants 

Maximum error at points 
within sampling instants 

Error in 

Error in 

Error in X^ 

Error in X^ 

Error 

Time of 
Occurence 

Error 

Time of 
Occurence 

Error 

Time of 
Occurence 

Error 

Time of 
Occurence 

Partial Matching 

H = [0 1], 

T * 2 sec 

'6.12X10 -3 

t = 7T 

-3.28xl0~ 4 

B 

6.12X10" 3 

15.5 sec. 

9.35xl0~ 4 

1 sec. 

Partial Matching 

H = [1 0], 

T = 2 sec 

1. 41*10~ 2 

t = 9T 

-1.66X10" 3 

t = 3T 

-1.41X10' 3 

17.5 sec. 

-1 . 66xl0~ 3 

6. 

Method of High 
Order Hold, 

T = 2 sec. 

N = 2 

5.3xl0~ 3 

t = 9T 

5x10 

t = 4T 

5.3xWT 3 

18. sec. 

-4 

5x10 

6,5 

Method of Multiple 
Sampling Periods 

T = 1. sec. , 

N = 2 

1, 267xl0 _3 

t = 9T 

1.308xl0~ 4 

t = 4T 

1.269X10' 3 

18.5 sec. 

-2.92xl0 -4 

J 

1. 
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4. DIGITAL REDESIGN OF THE CONTROL 
OF THE WOBBLE DYNAMICS OF THE SPINNING SKYLAB 

4 . 1 Introduction 

This report presents the results of the digital redesign of the control of 
a spinning Skylab space station. The dynamics of the spinning Skylab are first 
described. Only the wobble dynamics are considered however. The optimal 
feedback gains and the forward gains of the digital control are determined by 
use of the following point-by-point methods: 

a. Partial matching of states 

b. Exact matching at multiple sampling periods 

c. Exact matching with higher-order holds 

The simulation of the spinning Skylab system was carried out using the 
partial matching method and the multiple sampling period method. In practice, 
since the digital control as derived from the partial matching method is the 
easiest one to implement, and the simulation results are quite good, it is the 


recommended solution. 
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4 *2 Modeling of the Spinning Skylab 

A detailed description of the spinning Skylab space station can be found 
in the literature [ 3,4 ]. A simplified model for the spinning Skylab vehicle 
is shown in Figure 4-1 The variables of the system are indicated on the diagram. 

The mission of a spinning Skylab makes it necessary to point the 3-axis 
at the sun rather than to passively stabilize the steady-state rotation of the 
vehicle about its 3-axis. This will place the solar panels, which are lying 
in the 1-2 plane of the vehicle, normal to the impinging rays of the sun, 
making maximum use of solar energy. The control torques may be provided by 
three control moment gyros (CMg's). 

The dynamic equations of the system are expressed in vector form: 


•«. •• 


Mz_ + Dz^ + Kz_ = - v 

where the dots represent derivatives with respect to x = fit. 


(4-1) 


z. = W 1 <P 2 <P 3 m 1 u 2 y 3 ]’ 


(4-2) 


The control vector is given by 


V - v 2 v 3 v 4 v 5 v 6 r 


[ 


2 2 
Ijfl I 2 « 


1 / 


0 


0 


0]' 


(4-3) 
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- Y 3 

0 

0 

e 


l+K ] 

1-K„ 


0 

-e 

0 


1 

1 

0 


-SY, 


Y : 

1 

0 


5y, 


-V 1 


0 

0 

1 


Y 


(4-4) 


o -■ 

ll+K^ 

0 

0 

-2? 

0 


-(i+Ki> 


0 

0 

0 

0 

- 2 £ 


0 

0 

0 

0 

0 

2 


2 5 Yj 

0 

0 

0 

A, 


0 

2CY 3 

0 

-2Y, 

-2 

A, 


0 ) 

0 


V 3 


0 

0 

0 


(4-5) 


- K 1 ° 
1+K 

0 K 2 ( 'l-K 2 ' ) 


-tv. 


"Yi 


0 

0 

-5 


0 

0 

K 

o 


0 

0 

0 

0 


0 

0 


0 

0 


a 2 2 -l 


Yi(°3 +1) 

0 

0 

0 


(4-6) 
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If the value of is neglected, and £ = ^3 /|, ^2 = t * le mot i° n i s 

separated from the spin velocity motion. In other words, the differential equations 
represented by Eq. (4-1) are now divided into two sets of uncoupled equations. 

Wobble Dynamics 


*i - Vs ‘ (1 + K iV ' K A • V3 ■ ' v i 


1 +K 1 


f 1+K i 

1 

i i -k 2 J 

$ 2 + (1 + K 1 )^ 1 + K 2 

Lk 

( 1_k 2 J 


*2 - - v 2 


(4-7) 


(4-8) 


-yj>i + y 1 y 3 ~ y 1 <t > 1 + Y 1 <a 3 + Dy 3 = -v 3 

Spin Velocity Model 

*3 + Vl ‘ 2Y 3^2 ' - v 4 

»3 + S 1 + Vl ‘ 2 ^2 + a i 1J l " _V 5 


y* 2 + 2cf> 3 + 2y 1 + A 2 y 2 + (a 2 - l)y 2 = -v g 


(4-9) 


(4-10 ) 
(4-11) 
(4-12) 


In this study only the wobble dynamics will be considered. 

In order to obtain the state equations of the wobble dynamics, a state 


diagram is constructed in Figure 4-2 using Eqs. (4-7) f ( 4-8) * and (4-9) 
Applying Mason's gain formula to the state diagram of Figure 4-2 we have 




where it has been assumed that V 2 = 

The physical characteristics of tl 
1^ = 1.25 x 10 ^ kg-m^ 

I 2 = 6.90 x 10 ^ kg-m^ 

I ^ = 7.10 x 10^ kg-m^ 

r 1 = 0' 

^ = 23.3 m 
r 3 = -1.53 m 


m = 227 kg 
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spinning skylab are given below: 



k^ = = 146 N/m 

k 2 = 7.4 x 10 4 N/m 

d = d 3 = 0.04 (k 3 m) l/2 * 
1/2 

d 2 = 0.04 (k 2 m) ' = 164 

a = o.6 s " 1 




- 0.8478 


= 0.19717 


0.03471 


0.0534 


1.20 


A = 0.0534 

- 1.787 
• = 905.53 ’ 


1.787 


7.28 


Substituting the values of these parameters into Eq.(4-13) , we have 
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4 . 3 Optimal Control Laws for the Continuous -Pat a Controllers 

A total of three control laws were provided by NASA which represent the 
solution of the control of the wobble dynamics of the spinning Skylab vehicle 
by continuous-data controllers. These control laws are described as follows: 
Control Law No. 1. 


where 


T = a<}> + $to 


(4-15) 


1 - [T 1 

4 = [*! 

w = [ 


a = 


3 = 


to = 


T T T 
2 3 J 

K 0]* 


w r 

w 2 

“ll 

a l2 

a 21 

a a 

0 

0 

*11 

6 12 

t— l 
CM 

ca 

^22 

0 

0 

■1 

0 

0 

-1 

0 

0 


w 3 + fi] ' 


0 

0 

oj 


dj) 

dt 


0 n 0 

-R 0 0 

0 0 0 


4 


(4-16) 

(4-17) 

(4-18) 


(4-19) 


(4-20) 


(4-21) 


a 12 = e 


Bn - 1^6 


all other a . . = 0 
ij 

all other 3^ = 0 


Then 


T 1 ' “12*2 + B n”l = ( “l2 + B 11 B) *2 ' 


(4-22) 



(a) Linear case A: <S = - 1 . 27 , e = 5.9 

Then 

a ^2 = = 2.655 x 10 ^- 


3 n = 1^6 = - 0.9525 x 10 


Thus 


T 1 = 2.0835 x 10 6 £ 2 . + 0.5715 x lO 6 ^ 


v.. = j = 4 . 63 <j>„ + 1 . 27 <L 

I fi 

In terms of state feedback, = -G( 0 )j<, 


V;L - [0 4.63 0 1.27 0 0 ] 


y 3 

V ^ J 


(b) Linear case B • <S 23 - 1 , e = 2 . 
a = I l^ 2e = 0,9 X 106 


3 n = I fl 6 = - 0.75 x 10 


Thus , 

T = 0.45 x 10 6 i + 0.45 x lO 6 ^ 


V 1 = 




*1 + *2 


( 4 - 23 ) 

( 4 - 24 ) 


( 4 - 25 ) 


( 4 - 26 ) 


( 4 - 27 ) 


( 4 - 28 ) 
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v 1 o [0 1 0 1 0 0] 


(a 1 

h 

*2 

y 3 

♦l 

*2 

N 


Control Law No. 2 (Linear Quadratic Law) 


(4-29) 


The optimal control as obtained by the linear quadratic law is 


r , \ 
*1 


v 1 = [1.263 -0.777 -0.172 2.768 -1.792 -0.095] 


y 3 

v. 


(4-30) 



4. 4 Digital Redesign and Simulation 

The digital redesign of the control for thw wobble dynamics of the 
spinning Sky lab vehicle is performed by three different Poifjt-by-point methods : 

a. Partial matching of states 

b. Exact matching at multiple sampling periods 

c. Exact matching at higher order holds. 

In each case the feedback and forward gains are determined for the various 
control laws. The imulation is performed using the method of partial 
matching and the method of multiple sampling periods. 

From Eq. (4-14) the state equations of the system are 


x - Ax + Bu 


with x = c|> 2 y 3 


u = V;L 


^ 4> 2 v 3 3' 


(4-31) 


0 

0 

o 1 

1 

0 

0 1 

0 

0 

1 

0 | 

0 

1 

0 

0 

0 

o 1 

0 

0 

1 

0.0463 

0 

-0.4388 ( 
1 

0 . 

1.0462 

-,01313 

0 


I 

0 | 

-.1522 

0 

0 

1.0462 

0 

-3.22 | 

0 

1.0462 

-.06657 

{0 0 

0 

-1,245 

0 

-1.245] ’ 



The control vector u for the continuous system is 

u = E r - G x (4-32) 

— c c c— c 

where the feedback and forward gains of the continuous system for the three 


control laws are 
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Control Law 1A 


G = [0 
c 

4.63 0 1.27 0 

E = 1. 
c 


Control Law IB 


G = [0 
c 

1. 0 1. 0 0] 

E = 1. 
c 


Control Law 2 



G - [1.263 -0.777 -0.172 2.768 -1.792 -0.095] 

c 

E = 1. 
c 

The results obtained by performing the digital redesign of this system with the 
three different methods are described below, 
a. Partial Matching of States. 

The method of partial matching was used on control law 2 with a 

sampling period of T = 1 and a weighting matrix 

H = [1 1 1 1 1 1] 

The sampled- data gains are: 

G (1) - -.286054 T - 1 

s 

G (2) = .0711766 

s 

G (3) = .161217 

s 

G (4) = -.946105 
s 

•G g (5) = -.113347 

G (6) = .121896 

s 

E (1) = .246084 

s , 

The continuous and sampled-data systems were simulated and the results 
are shown in Figure 4-3 '■ through 4-9 ; The states and controls, and 

the errors between the continuous and sampled states and controls 
are shown. The results show that for T = 1 the point-by-point method of 
partial matching yields an acceptable redesign. 



Time 


Feedback Gains from States 


Forward Gain 


Interval 

X 1 

X 2 

X 3 

X 4 

x 5 

X 6 

r 

.0 

-4.92398 

-801.398 

55.8475 

-826.826 

. -283.221 

9.22254 

830.723 

.5 

-14.5192 

. -594.385 

, 153.307 

-549.619 

-706.054 

-1.40256 

546.516 

1.0 

-13.3919 

-867.241 

17.1915 

-1058.39 

37.3406 

74.8695 

870.569 

1.5 

40.4243 

-619.488 

68.0836 

-342.295 

-929.208 

-51.2130 

291.012 

2.0 

35.3183 

-859.229 

85.9373 

-1008.44 

-64.9085 

22.6930 

597.738 

2.5 

.156333 

-624.771 

83.1258 

-662.564 

-514.606 

22.8996 

516.595 


-3 -3 

Table 4-1 Forward gain E X10 and Feedback gains G X10 for the sixth 

s s 

order model of the wobble dynamics of the Spinning Sky lab. 
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b. Exact Matching by Multiple Sampling Periods. 

The method of Multiple Sampling Periods was used on control law IB 

for T = 0.5. That is the gains are changed every 0.5 seconds and all 

the states are to be matched every 6 x 0.6 = 3 seconds. The sampled- •• 

data gains are given in Table 4-1 for G and E . The continuous and 

s s 

sampled systems were simulated and the states, controls, and the errors 
between the states and controls of the sampled and continuous systems 
are shown in Figures 4-10 to 4-16. The states of the two systems are 
reasonably well matched considering the size of the sampling period. 

c. Exact Matching with Higher-Order Holds. 

This method for digital redesign has been described in Chapter 3. 

With n/m = 6 , it appears that a 5th-order hold is adequate. Thus, 

N = 6 is used and the coefficients of the various terms of the gain 
matrices are calculated with T = 0.2. Table 4-2 lists the gains for control 
law 1A. 

The simulation has not been performed with this method of digital 
redesign, but it is expected that the results will be similar to those 
obtained by the method of multiple sampling periods. 



Table 4-2 


Coefficients 


of 

X 

1 

*2 

1 

.29116 

-464.55 

t - kT 

-.57297 

10.953 ' 

RSHESSK 

IH 

.41414 

-3,52 

(t - kT) 3 
3! 

-.18961 

2.1888 

(t - kT) 4 
4! 

.48 

-6.922 

(t - kT) 5 
5! 

-.5382 

9.699 


Feedback and Forward Gains for the 
wobble dynamics of the Spinning Skylab , 

Control Law 1A. Method of Higher-Order Holds, 


Feedback Gains : 
Feedback from States 


X 3 

X 4 

X 5 

-.2666 

-127.72 

2.3438 

'1.0845 

3.4844 

-10.7695 

-,45875 

-1.36 

4.29 

.2446 

.7616 

-2.3768 

-.68544 

-2.2682 

6,8147 

.88064 

3.0556 

-8.966 


Forward Scale 

Gains Factor 

X 6 

-.08154 100.56 X10“ 2 

.17682 -2.75 

-.13752 1.05 X10 2 

.07324 -.5888 X10 4 

-.2046 1.6947 X10 5 

.26099 -2.2364 X10 6 


















































3. oo 40.00 80.00 ' 120.00 160.00 200.00 240.00 280.00 

T(X1C 

Figure 4-3 State trajectories x.^(t) and x.^(t) for the 6th order model of the 

Spinning Skylab Satellite. Digital redesign by the point by point 

method of partial matching. 
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Figure 4-4 


80.00 120.00 . 160.00 200.00 240.00 280.00 

tCxicT 1 ) 

State trajectories x 2c (t) and x (t) for the 6th order model 
of the Spinning Skylab Satellite. Digital redesign by the 
point by point method of partial matching. 




UO.OO 80.00 120.00 160.00 200.00 240.00 

Figure 4.5 State trajectories x^Ct) and x^ g (t) for the 6 th order model 
of the Spinning Slcylab Satellite. Digital redesign by the 
point by point method of partial matching. 


280.00 

T(X10 -1 ) 




Q-QO 40.00 80.00 120.00 160.00 200.00 240.00 280.00 

T (X10 -3 ) 

Fig. 4-7 State trajectories x (t) and x c (t) for the 6th order model of the 

->c ->s 

Spinning Skylab Satellite. Digital redesign by the point by point 
method of partial matching. 


200.00 


240.00 


Figure 4-8 State trajectories x^t) and x (t) for the 6th order model of the Spinning Sky lab 
Satellite. Digital redesign by the point by point method of partial matching. 
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, . .i ' ?—». ' ^ 

Figure 4-9. Control trajectories u (t) and u (t) for the 6th order model 

C S 

of the Spinning Skylab Satellite. Digital redesign by the i 
point by point method of partial matchint. 
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Figure 4-10 State trajectories x^Ct) and*-x (t) for the 6th order model of the Spinning Sky lab 

Satellite. Digital redesign by the method of exact matching at multiple 
sampling periods. 
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Figure 4-12 State trajectories x„ (t) and x 0 (t) for the 6th order model of 

jc 3s ' / 

the Spinning Skylab Satellite. Digital redesign by the method of 
exact matching at multiple sampling periods. 


3, gg 40.00 80.00 120.00 180.00 200.00 240.00 

Figure 4-13 State trajectories x. (t) and x. (t) for the 6th order model of: 

4c 4s 

the Spinning Skylab Satellite. Digital redesign by the method 
of exact matching at multip^^ sampling periods. 
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Figure 4-14. State trajectories x c (t) and x c (t) for the 6th order model 
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of the Spinning Skylab Satellite. Digital redesign by the 
method of exact matching at multiple sampling periods. 
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Figure 4-15 State trajectories x, (t) and x.. (t) for the 6th order 
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model of the Spinning Skylab Satellite. Digital redesign 
by the method of exact matching at multiple sampling periods. 
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Figure 4-16 Control trajectories u (t) and u (t) for the 6th order 

c s 

model of the Spinning Slcylab Satellite. Digital redesign 1 

by the method of exact matching at multiple sampling periods. 
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5. OPTIMAL REGULATION OF THE TWELFTH-ORDER SPINNING SKYLAB 
5 . 1 Introduction 

This report considers the calculation of the optimal feedback 
gains and simulations of the closed- loop system of the twelfth-order 
spinning skylab . The linear quadratic law (LQL) method is used for 
the optimal design of the skylab system. 

It has been determined that the overall skylab system is 
uncontrollable with only one input to one of the wobble, axes. An 
additional input to one of the spin axes has been found necessary 
to provide controllability of the system. 

The optimal feedback gains have been obtained by feeding back 
eleven of the states, with ti<ro different weighting matrices for Q. 

Simulation results are presented for the twelfth-order system 
using the optimal gains as obtained from the design of the sixth- 
order wobble dynamics, as well as using the optimal gains designed 
with the eleventh-order model. The eleventh-order model was used by 
deleting the state which is the cummulative deviation of the 
position of the spin axis. Since it is of interest only to regulate 
9 ^ about its zero reference, the state <p^ is not necessary in the 
optimal design. 



122 


5.2 Modeling of the Twelfth-Order Spinning Skylab 

The dynamic equations of the spinning skylab space station are expressed 
in vector form: 


Mz + Dz + Kz^ = - v 

where the dots represent derivatives with respect to T 

* ■. 

£ ■ [* x <P 2 <P 3 v*! p 2 y 3 3 

The control vector is given by 


V =[„1 v 2 v 3 v 4 v 5 v 6 ] 


Sit 


(5-1) 


(5-2 ) 


T 


[■ 


3 2 

1 1 « i 2 n 


‘3 


IJT 


0] 


(5-3) 


0 


0 
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. shall consider the wobble motion together with the spin velocity motion 
so that Eq. (5-1) represents a twelfth-order system. 

In order to obtain the state equations of the overall system, a state 


diagram is 

constructed as shown in Figure ! ( 5 — 1 ) The state equations 

are 


determined 

from Figure( 5-l)by use of Mason's gain formula. We have 




*1 


O 

« 

• 

• 

• 

0 

t — i 

0 

• 

• 

• 

• 

0 

0 

\ 


r 0 


*2 

. 

0 0 . . . 0 i 0 1 . . . . 0 

*2 


0 


*3 


■ ' - I * 

*3 


0 


*1 


. . 1 

'■i 


0 


^2 


0 0 • • 0 } 0 0 . . . 1 0 

V 


0 



- a 

0 0 . , . • 0 j 0 0 . , . . 1 

1^ 

P 3 


0 


h 


1 


+ 

! 

! 

- d l 


h 


1 

^2 

: 

■ 

0 


h 


l 

-M _1 K } -M _1 D 

*3 


0 


h 



• 

^1 


0 


*2 


I « 

d 2 




•• 



♦ 



> 

y 3 



% 


- d l 



/ 

. 4 

[ 

* 





(5- 

-7) 





125 


where 


M 


-1 


M 1 K 


M _1 D 


d l 

0 

d l 

0 

0 

-«Vi ' 

0 

(1-Y 3 )cd 2 

0 

-^i cd 2 

5Yl d 2 C 

0 

0 

^3 ?cd 2 

0 (l-? 2 Yl c)d 2 

" Y 3 d 2 

0 

0 

5cd 2 

0 

" d 2 

d 2 

0 

^ d l 

0 

-Edi 

0 

0 

(1 -V d l 

/ 

1 

rH 

-d 

0 

' 52)d i 

0 

0 

-«Vi 

J 

r 

a 71 

a 72 

a 73 

a 74 

a 75 

a 76 

00 

H 

a 82 

a 83 

a 84 

a 85 

a 86 

a 91 

a 92 

a 93 

a 94 

^5 

^6 

a 10,l 

a 10 ,2 

a 10,3 

a 10 ,4 

a 10,5 

a 10,6 

a ll,l 

a ll,2 

a ll,3 

a ll ,4 

a ll,5 

a ll,6 

a 12,l 

a 12,2 

a 12,3 

a 12 ,4 

a 12 ,5 

a 12 ,6 

*77 

a 78 

a 79 

a 7 ,10 

a 7,ll 

a 7,12 

a 87 

a 88 

; a 89 

a 8,10 

a 8,ll 

a 8 ,12 

a 97 

a 98 

a 99 

a 9,10 

a 9,ll 

a 9 ,12 

a 10 ,7 

a 10 ,8 

a 10 ,9 

a 10 ,10 

a 10,ll 

a 10,l 

a ll,7 

a ll,8 

a ll,9 

a ll,10- 

7 a ll,ll 

a ll,l 

a 12 , 7 

a 12 ,8 

a 12 ,9 

a 12 ,10 

a 12 ,11 

a 12,l 


( 5 - 8 ) 


( 5 - 9 ) 


( 5 - 10 ) 



126 


1 1 - Yl (l + S 2 ) 
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5 . 3 Controllability of the Skylab System 
Given a linear time- invariant system 

x = Ax + Bji (5-14 ) 

where x is n-dimensiona.1, and is m-dimensional. A common method 
of checking the controllability of the system is to see if the 
following rank condition is satisfied; 

Rank [B AB A 2 B ... A n " 1 B] = n ' 

. (5-15) 

2 

For the skylab system with one input, » the above 

matrix is square and thus the rank condition can be checked by 
calculating the determinant of the matrix. Using the numerical values 
of the system parameters given in Chapter 4, section 4.2 i 

. * " I ‘ 

of this report, Ithe determinant of the matrix of 

1 8 

Eq. (5-15) is found to be 4.46 x 10 . Therefore, it would appear 

that the system is completely controllable. _ However, due to the 
numerical difficulties encountered in calculating the Riccati gains 
as well as in shifting the closed-loop eigenvalues, it was decided 
to perform an alternate check on the controllability of the system. 

An alternate test on the controllability can be made by performing 
the similarity transformation 

x - P£ ( 5-16) 

on Eq.(5-14) where P is an-n x n nonsingular matrix. The new system 
of state equations are written 

| (5-17) 


v = A 2 . + r J± 



where 


132 


A = P -1 AP 

and 

r = p _1 b 


(5-18) 

(5-19) 


For A with r pairs of distinct complex eigenvalues and n - 2r 
distinct real eigenvalues, A is written 



with 0 ± ± jug, i = 1, 2, ... r, denoting the ith-pair of complex 
eigenvalues; and , i = 2r+l, ... . , n denotes the n-2r real 
eigenvalues . 

For complete controllability, the rows of T which correspond 
to the distinct real eigenvalues must not contain all zeros, and at 
least one of the two rows which corresponds to the rows of S_^ for 
i = 1, 2, ..., r, does not contain all zeros. 

For the eleventh-order system, the eigenvalues of A are given 
in Table 5-1 . 

The T matrix which corresponds to the same order as these 
eigenvalues is obtained as 



133 


* 

-8.16 

X 

io -2 

1.67 

X 

10“ 3 

-9.32 

X 

io" 1 

2.99 

X 

io -2 

-1.32 

X 

IO" 3 

9..0 

X 

io" 2 

-1.76 

X 

io" 1 

1.54 

X 

io- 10 

1.97 

X 

io' 3 


-3.19 

-2.79 x lO -16 ^ 

Since the last element of T is practically zero, the system 
is uncontrollable. In fact, the state which corresponds to the 
eigenvalue at the origin, s = 0 , is the uncontrollable one. This 
explains the difficulty in attempting to move this eigenvalue into 
the left-half plane. This alternate controllability test also shows 
the interesting fact that in practice, when there are very large 
and very small parameter values in a large-scale system, the control- 
lability check using Eq. ( 5-15 ) may be unreliable, due to the repeated 

products of very large and very small numbers . 

2 

When the control V. = T„/I„fi is used in addition to the control 

4 3 3 

V^, carrying out the similarity transformation, we have 
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-8.2 

X 

io“ 2 

-1.2 x 
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-O* 
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io“ 3 

-7 x 

t-O 
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»H 

-9.3 
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rH 
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• 

ro 
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io -3 

3 
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io -2 
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10“ 2 

-1.3 

X 

IQ ' 3 

1.03 


9.0 

X 

io -2 

-2 x 

10 -2 

-1.8 

X 

io" 1 

-2.9 x 

io- 10 

1.5 

X 

io“ 10 

3.6 x 

K.r 9 

1.9 

X 

io -3 

1.6 x 

io* 2 

-3.2 



1.2 x 

io' 4 

-2.8 

X 

io -16 

9.9 x 

io -1 


(5-23) 


Now since all the elements of the last row of T are not zero, 
and not all rows which correspond to the complex eigenvalues are 
zero, the system is completely controllable. 
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5.4 The Eigenvector Method of Calculating Riccati Gains 

Consider the matrix quadratic equation. 

- KA - A'K + KSK - Q = 0 ; (5-24) 

where the coefficient matrices A, S, and Q are given n x n matrices. The 
solution K is also an n x n matrix. 

It will be shown later that(5-24)possesses many solutions one of which is 
positive definite, one negative definite and the rest indefinite. We will be 
concerned with the positive definite solution which is denoted by K + . 

Let us form the 2n x 2n matrix 



(5-25) 


and denote the 2n eigenvalues of M by 1,, i = 1, .... 2n and the 2n eigenvectors 
of M by a_^, i = 1.... 2n. Partitioning the eigenvectors we have 


r b 


(5-26) 


where a is a 2n vector, b is an n vector and c is an n vector, 
i i i 

The following properties which are proven in [5J form the basia for the 
computation of K. 

i) If is an eigenvalue of M, there is an eigenvalue Xj of M such that 

* 

Aj = i»j = i>2, . . . ,2n (5-27) 

where the asterik denotes complex conjugate. 

ii) M has no purely imaginary eigenvalues. 

iii) If the coefficient matrices A, S and Q are real, eigenvalues 

off of the real axis appear in conjugate pairs 
* 

■ \ m *i 
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The eigenvectors corresponding to the above eigenvalues 
also appear in conjugate pairs. 



iv) The n eigenvalues of M with negative real parts correspond to 

+ 

the eigenvalues of (A - SK ) . 

If a^, , ...» a n are the eigenvectors corresponding to the n 

eigenvalues with negative real parts , then 

K + = [b 1 b 2 ... b n ] [c 1 c 2 ... c q ] 1 ; (5-28) 

Every solution of F(K) » 0 is the form 

R+ [b 1 b 2 ... b^]^ c 2 ... cj" 1 ; (5-29) 

If the inverse exists, where the eigenvectors correspond to any 
n eigenvalues of M. 

Any linear combinations of the vectors b.. , b_, ..., b can be used 

1 Z n 

in(5-28) and ( 5-29 )without affecting the result if the same 

linear transformation is applied to the vectors c. , c. , ..., c . 

12 n 

The computer algorithm then consists of the following steps: 

i) Form the M matrix from A, S, and Q. 

ii) Calculate the eigenvalues and eigenvectors of M. A University of 
Illinois subroutine using the QR iteration technique is used 
for this calculation. > 

iii) Using a sifting routine separate out the n eigenvectors of M 

corresponding to the eigenvalues with negative real parts. Form 

- ... . ..j* 

a real 2n x n matrix, D, from these eigenvectors. For complex 
pairs one column of D consists of the real part of the eigen- 
vector and one column consists of the imaginary part. 


v) 


vi) 


vii) 
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iv) Separate the D matrix into two n x n matrices, B and C, as 
shown in (5-26) 

v) The solution is given as 

K + = B(f 1 (5-30) 

Note that the eigenvectors a^ either appear as purely real or in conjugate 
pairs. Therefore it is always possible to take linear combinations of a^, to 
form the real matrix D. To see this let E be the 2n x n matrix consisting 
of the n eigenvectors of M corresponding to eigenvalues with negative real 
parts. Let the n x n matrix G be the complex linear transformation on E which 
results in the matrix D. 



For example let the first column of E be a real eigenvector and the next 

two columns consist of a complex pair of eigenvectors. The G matrix would then 

be of the form 



If we partition E into two n x n matrices E^ and E 2 it can be seen that the 

+ 

final solution, K , does not depend on the transformation G. 



K + = BC _1 - E 1 G(E 2 G)~ 1 = EjE" 1 


- • t * 



Numerical Results 
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5.5 


The eigenvalues of the uncontrolled eleventh order open loop 
system are shown in Table; 5-1 Note here that in all cases the 
coupling parameter between the spin and wobble dynamics in nonzero 
I* = 1.56. By investigating the eigenvectors associated with these 


eigenvalues, it can be seen that major contributions of the zero 
eigenvalue is to the state ^ . Therefore the use of a control law 
which cannot shift this zero eigenvalue will result in a steady 
state error in (*>„ 

J • 


The eigenvalues of the eleventh order closed loop system using 
control law 2 given in Chapter 4, Eq. (4-16) 

are shown in Table 5-1. Since this is a linear quadratic control law 
for the regulation of the wobble dynamics of the Spinning Skylab , it 
might be expected that it will not adequately control the spin 
dynamics. This is indicated in Table 5-1 by the continued presence 
of the zero eigenvalue. By examining the eigenvectors corresponding 
to the closed loop eigenvalues we may again see that the major 
contribution of the zero eigenvalue is to the state 4>^. 

Since it is shown in the last section that the zero eigenvalue 
mode is uncontrollable from the input , the additional input 
is used when calculating the optimal Riccati gains for the eleventh 
order system. The Riccati gains are calculated for the following 
two cases of weighting matrices: 

Case 1 R = I, Q = I 

Case 2 R = I, 0 = diagonal [10 10 1 1 1 10 


10 1 1 1 1 ] 
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The feedback gains for these two cases are shown in Table 5-2 The 
closed loop eigenvalues for case 1 and case 2 are shown in Table 5-1 
In both cases the zero eigenvalue has been shifted into the left half 
plane resulting in an asymptotically stable eleventh order system. 



THE EIGENVALUES OF THE OPEN AND CLOSED LOOP SYSTEMS 


Open Loop 
A Matrix 

Control Law 2 
(LQL design of 
wobble dynamics) 

Case 1 

Case 2 

-0.603 ± j 30 . 15 

-.604 ± j 30. 15 

-0.604 ± j 30 . 15 

-0.604 ± j 30 . 15 

-.0328 ± jl.76 

-.048 ± j 1.68 

-0.369 ± jl.56 

-0.159 ± j 1.64 

-.0276 ± jl.36 

-.0277 ± jl.36 

-0.52 ± jl.24 

-0.516 ± jl.24 

0.0 ± jl.O 

-.078 ± j0.92 

-0.080 ± j0.938 

-0.077 ± j0.92 

-.00056 ± j0.29 

-2.53 

-0.974 ± j0.653 

-3.82 

0.0 

-0.611 

-1.01 

-1.06 


0.0 


-1.01 


Table 5-1 The eigenvalues of the open and closed loop 
for the eleventh order skylab. 


0<7l 











States 

CASE 1 

CASE 2 

Control 

Control V, 
4 

Control 

Control V. 

4 

1 

-.0513 

-.043 

-2.67 

-.0054 

2 

.030 

-.053 

2.15 

-.09 

3 

-.009 

.930 

.041 

.9 30 

4 

-.094 

-1.82 

.077 

-1.83 

5 

-1.29 

.052 

-1.28 

.0254 

6 

-1.78 

-.016 

-4.4 

-.0176 

7 

.828 

-.058 

3.27 

.175 

8 

.023 

-1.0 

-.001 

-1.0 

9 

.034 

.918 

.015 

0.919 

10 

.0021 

.0043 

.0338 

.0043 

11 

-.46 

.0082 

.155 

.01 


Table 5-2 The optimal Feedback Gains for 

the eleventh order Skylab. 
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5 . 6 Simulation 

The 12th order model for the Spinning Skylab is simulated with 
the following two feedback gains 
a. 


G 

c 


-.051 

.03 0 

-.009 -.094 

-1.29 -1.78 .828 

.023 

.034 

.002 

-.043 

-.053 0 

.930 -1.82 

.052 -.016 -.058 

- 1.0 

.918 

.004 


This gain is obtained by solving the Riccati equation for the 
eleventh order system with two control inputs and with the weighting 
matrices Q = [I] and R = [I], (see Case 1 in Table 5-2 ). 


G = [-1.263 .777 000 .172 -2.78 1.792 000 .095] 

c 

This gain is obtained by applying the linear quadratic law (LQL) 
method to the 6 th order wobble dynamics of the Spinning Skylab, 

In both simulations the initial states are 


x( 0 ) = [.01 .01 0 0 0 0 0 0 .01 0 0 0 ]’ 

Figures 5-2 through 5-9 show the state and control trajectories for 
case a and Figures 5-10 through |5-17show the state trajectories for 
case b. 

As can be seen from the simulation results, in case a all the 
states are regulated about the zero reference except which undergoes 
a constant shift in its steady state magnitude. In case b all states, 
except x 3 and x^ , are regulated about the zero reference; x 9 undergoes 
a constant shift in its steady state magnitude and x^ continually 


-.46 

.008 


increases. 



State x (t) (X10 


CD 

CD 



Figure 5-2 State trajectory x^Ct) of the 12th order 
model for the Spinning Skylab Satellite. 


T(xlO _1 ) 


143 


State x„(t) (X10 


-.051 .03 


-.009 -.094 -1.29 -1.78 .828 .023 .034 .002 -.46 


-.043 -.053 0 .93 -1.82 .052 


6 -.058 -1. .918 .004 .008 


x(0) = [ .01 .01 0 0 0 0 0 


0 0 ]' 
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Figure 5-3 State trajectory (t) of the 12th order 
model for the Spinning Skylab Satellite. 
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State trajectory x^(t) of the 12th order model for 


Spinning Skylab Satellite 
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Figure 5-5 State trajectory x^(t) of the 12th order 


model for the Spinning Skylab Satellite. 
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State x_(t^ (xlO 
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Figure 5-12 State trajectory x^(t) of the 12th order 
model for the Spinning Skylab Satellite. 
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Figure 5 16 Control trajectory u(t) for the 12th order model of the Spinning Sky lab. Feedback 
gain obtained from the LQL design of the wobble dynamics (6th order model). 
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DIGITAL REDESIGN OF THE 11th ORDER MODEL OF THE SPINNING SKYLAB 


The design of the continuous system for the 11th order model of the 
Spinning Skylab was performed by use of the LQL method. The feedback 
gain matrix for Q = I, R = I was found to be (see Chapter 5, Table 


5-2). 



1 

1 







-.513 

.302 

-.009 

-.094 -1.29 -1.783 

.828 

.0229 

.034 

.0021 

-.457| 

c 

-.043 

-.053 

.93 

-1.82 -.0515 -.016 

-.058 

-1.0 

.918 

.0043 

. 0082 J 


The digital redesign of this system is performed by use of the point by 
point method of partial matching. The weighting matrix H is chosen as 


H » 


f 

1 

0 


1 

0 


0 

1 


0 

1 


1 1 
0 0 


1 

0 


0 0 0 1 
1110 . 


and the feedback gains for the digital system are determined for two 
different sampling periods, T = .2 and T = .8. The choice of H was 
made such that matching is obtained for the sum of the states of the 
wobble dynamics as well as the sum of the states of the spin dynamics 
for each sampling instant ( See Chapter 1, Section 1. ) 
for a discussion on how H should be chosen). With T = .2 the feedback 
gain matrix of the sampled data system is 


G = 

r -.442 

.183 

-.011 

-.0495 

-.777 

-1.41 

.450 

.0178 

.019 

.031 

-.434' 

w 

with 

(-.043 
T = .8 

-.025 .696 

we have 

-1.53 

-.013 

-.063 

-.014 

-.875 

.85 

.015 

.004. 

G = 

- . 332 

.095 

-.0045 

.0055 

.128 

.775 

.21 

.006 -. 

0057 

.0017 

-.22 

w 

-.032 

.026 

-.0015 

- . 66 

-.125 - 

.137 

066 

-.51 

548 

-.009 

-.025 , 


The simulation of the 12th order model of the Spinning Skylab is performed 
with each of the above gains. Figures 6-1 through i.6-14 show the results 
with T = 0.2. These figures show the state, control and errors in state 



and control trajectories. Figures 6-15 ! through 6-22; show the 
results for T = 0.8. 

It is apparent that the point by point state comparison method of 
partial matching yields acceptable redesign with T =» .2 as well as T = 



0.00 40.00 80.00 120.00 160.00 200.00 240.00 280.00 T (X10 ) 

Figure 6-1 . State trajectories x^t) and x (t) for the 12th order model of the Spinning Skylab 
Satellite. Digital redesign by the point by point method of partial matching. 




80.00 120.00 160.00 200.00 240.00 280.00 

T(X10 -1 ) 

Error in state trajectories, x^Ct) - x (t) , for the 12th order 
model of the Spinning Sky lab Satellite. Digital redesign by the 
point by point method of partial matching. 
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• * _ "1 

Figure : 6-3 State trajectories ^(t) and x 2s (t) for the 12th order x(X10 ) 

model of the Spinning Skylab Satellite. Digital Redesign 

by the point by point method of partial matching. 
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Figure ' 6-4 Error in state trajectories, x^t) - x^Ct), for the 12th order 

model of the Spinning Skylab Satellite. Digital redesign by the point 
by point method of partial matching. 
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Figure 6- 
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120.00 160.00 200.00 2 * 40.00 

State trajectories x^Ct) and x^ (t) for the 12th order model 
of the Spinning Sky lab Satellite. Digital Redesign by the point 
by point method of partial matching. 
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6-9 State trajectories x (t) and x. (t) for the 12th order model of 

jc 5 s 

the Spinning Skylab Satellite. Digital Redesign by the point by 
point method of partial matching. 


Error in states, x. (t) - x_ (t) , (X10 
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Figure, 6-10 Error in state trajectories, ^(t) - x 5g (t), for the 12th order 
model of the Spinning Skylab Satellite. Digital redesign by the 
point by point method of partial matching. 
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QO UQ.OQ 80. oa 120.00 160.00 200.00 240.00 280.00 

! '• 

Figure 6-12 . Error in state trajectories, * 6c (t) - x 6s (t), for the 12th order T(X10 ) 

model of the Spinning Skylab Satellite. Digital redesign by the 
point by point method of partial matching. 
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Figure 6-13 Control trajectories u lc (t) and u^Ct) for the 12th order model 
of the Spinning Skylab Satellite. Digital Redesign by the point 
by point method of partial matching. 
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0.00 40.00 80.00 120.00 160.00 200.00 240.00 280.00 

Figure 6-15 State trajectories ^(t) and x (t) for the 12th order model T(X10 _1 ) 

of the Spinning Skylab Satellite. Digital Redesign by the 
point by point method of partial matching. 
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Figure 6-16 State trajectories x^Ct) and x (t) for the 12th order model of the Spinning Skylab 
Satellite. Digital Redesign ty the point by point method of partial matching. 
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Spinning Skylab Satellite. Digital Redesign by the point 
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Figure 6-22 Control trajectories u. (t) and u„ (t) for the 12th order 

2 c 2s 

model of the Spinning Skylab Satellite. Digital Redesign by 
the point by point method of partial matching. 
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7. ANALYSIS OF THE ATTITUDE DYNAMICS OF A 
LARGE SPACE TELESCOPE (LST) MODEL 

7.1 j Introduction 

The problem considered in this report is the determination of the condition 
of limit cycle in the Large Space Telescope (LST) using Control Moment Gyros 
(CMG's). Due to the on-board digital control system, the dynamics of the LST 
are represented by the sampled-data system of Figure 7-1 j 

A suitable model of the CMG is considered to be a combination of a nonlinear 
element and a linear element [6] . j The simplified input-output relationship of 
the nonlinear element is shown in Figure 7-2 . . 

One established method of predicting limit cycles in nonlinear sampled- 
data systems is the discrete describing function method [7].j The stability 
study of the LST system with k = 0 for the nonlinearity in Figure 7-2 can be ! 
readily predicted [12]. j This report is concerned with the derivation 

of the discrete describing function when k 4 0 in Figure 7- 2; • 


below: 


Specific system models have been considered by NASA, and these are listed 

7 : 

Case I. (First-order CMG, perfect attitude rate sensor) 

H(s) ° 1 - perfect attitude rate sensor 


w 


G(s) - 


s+w 


first-order CMG 


Case II. (Second-order CMG, perfect attitude rate sensor) 
Ii(s) = 1 perfect attitude rate sensor 


K, 


G(s) 


(7-1) 


s + 8 1 S+ S 0 


(7-2) 




Figure 7-1 ‘ 


I = Vehicle moment of inertia about the chosen axis 
Kq = control gain for 0 

= control gain for 0 

H(s) = transfer function of attitude rate sensor f rate gyro 
G(s) = linear transfer function of CMG 
N = nonlinear portion of CMG 
X = command 

6 = actual vehicle attitudes 
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Case III. (First-order CMG, first-order attitude rate gyro) 


to 


H(s) 


s+to 


imperfect attitude rate sensor 


to 


G(s) - 


s+to 


Case IV. 


(Second-order CMG, first-order rate gyro) 


H(s) 


to 


s+to 


G(s) = K G 
2 

s +g 1 s+g Q 

List of System Parameters 

I = 1.356 x 10"* Nms 
K Q = 4.6843 x 10 6 Nm 

K x - 1.1153 x 10 6 Nms 

g 0 = k g = 1q4 (rad/s) 2 

g^ = 141.4 rad/s 

to =70.7 rad/s 
c 

to = 30 rad/ s 
g 


(7-3) 


(7-4) / 


(7-5) 

(7-6) 


For the purpose of stability analysis, the system shown in Figure 7-1 tan 
be represented by the equivalent system of Figure 7-3 ' The two systems are 
equivalent in the sense that the characteristic equations are identical. The 
equivalent transfer function of the system shown in Figure 7-3 is given by 


K +K H(s)s ' 

G e (s) = ^ G(s) (7-7) 

Is - 



Thus, for the four cases listed previously. 


Case I. 


G e (s) 


( K l s + K o)u c 

Is 2 (s+W ) 
c 


Case II. 


Case III. 


Case IV. 


G e (s) 


G (s) 
e 


G e (s) 


< k i sH V k g 

S3 — - - ...... 

2 2 

Is (s +g 1 s+g Q ) 

_ V+V/K.y 

2 W 

Is (s-k/J ) (s+ 0 ) ) c 
c g 

K o s+K o tJ g +K i U) g s 

2 2 

Is (s +g 1 s+g 0 )(s+w g ) 


186 


(7-8) 

(7-9) 


(7-10) 


(7-11) 


7.2 . The Discrete Describing Function (DDF) 

The discrete describing function may be used to determine the condition 

of sustained oscillation in a nonlinear sampled-data system of Figure 7-3 It 

is assumed that the reference input x(t) is zero, and the actuating signal e(t) 

is sinusoidal with period = nT, where n = 2, 3, 4, ..., and T is the sampling 

period in seconds. This assumption is justified, since, in practice, G (s) has 

e 

the characteristics of a low-pass filter. 

From the analytical standpoint, an equivalent system of Figure 7-4 is used 

; > 

to replace the block diagram of Figure 7-3 In Figure 7-4 the zero-order hold 
(zoh) and the nonlinearity are transposed. This does not affect the system 
behavior if the nonlinear element is amplitude dependent only. However, the 
nonlinearity must be redefined to act on the strength of impulses. 

The closed-loop transfer function of the system of Figure 7-4 is 


(7-12) 


6(z) N(z) G(z) 
X(z) " l+N(z)G(z) 









where 


G(z) 


z 


f 

1-e 


-Ts 


s 


G (s) 
e 


= (1 - z 1 )s[G e (s)/s] ■ (7-13) 

and N(z) is defined as the DDF; i.e., 

N (z) = v(z)/K(z) (7-14) 

As in the conventional describing function technique, self -sustained 
oscillation occurs in the system when the following equation is satisfied: 

1 + N(z)G(z) = 0 (7-15); 


or 

G(z) =-l/N(z) (7-16) 

The derivation of the DDF of a nonlinear element is based on the assumption 
that the input to N is a sinusoidally amplitude-modulated impulse train. The 
output of N is also a train of impulses, but with the strengths of the impulses 
determined by the nonlinear characteristics. 

Let the input to the sampler be 


e(t) = E cos (w t +6) 
c 

The z-transform of the last equation is 

Ez [ (z- cos to T) cos <f> - sin W T sin <J>] 
E(z) = 2 5 9 

z - 2z cos w T + 1 
c 


(7-17) 


_(7-18) 


The period of the self-sustained oscillations is designated by T , where 

T c = 27 t/oj c = NT (7-19); 

and N = 2, 3, 4, ... . In addition to the characterization of the oscillations 
by T c> we need to classify the response with respect to the waveform of h(t) 
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which is the output of the zoh. Since e(t) is a sinusoid, it can be shown that 


for all T c = NT, N = even integers, the number of positive relay corrections 
is equal to that of the negative relay corrections during one period T^. This 
number of correction is designated by A = (i, j), where i is the number of 
positive relay corrections, and j is the number of negative relay corrections, 
during each period T^. For T^ = 2T, it is apparent that A can only be (1,1). 
For T £ = 4T, A = (1,1), or (2,2). For N = odd integers, the number of 
positive and negative relay corrections may differ by one. Therefore, in the 


representation for A, i and j may not be equal. 


For a given period T^, the loci of - 


form regions in the magnitude 


versus phase domain, each for one possible A or A,^. These regions are 

defined as the critical regions for the specific nonlinearity, T , and A. The 

c 


symbols , - - 


and 


max 


are used to indicate the boundaries of the 


mm 


critical region. The condition of self-sustained oscillation in the system of 
Figure 7-4 is determined with the graphical procedure described as follows: 

1. Plot G(z) in the gain-phase (db-degrees) coordinates for 
T^ = 2T, 3T , 4T , ... using T as a parameter on the loci. 

2. Superpose on the G(z) plot the family of critical regions of 


- ^rrr for T = 2T, 3T , 4T, ... . 

N(z) c 

3. If the portion of the G(z) locus for some T falls within the 

c 

critical region of the same T^, then there exists a set of 
E, <P, A, such that Eq. (7-16) is satisfied. Consequently, 
self-sustained oscillations characterized by the A, e, , and T 
will occur for the range of T which corresponds to the portion of 
G(z) found inside the critical region. 
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4. The sampling periods T along the portion of G(z) locus outside the 
critical region of the same T^ and some A correspond to operations 
which will not give rise to self-sustained oscillations of the 
prescribed mode. 

The typical procedure of deriving the critical region for the nonlinearity 
of Figure 7-2 for k ^ 0 is described in the following section for T = 3T, or 
N = 3, for all possible combinations of A. The procedure for generating the 
gain-phase plots of G(z) is described in a later section. 


7.3 : Derivation of the Critical Regions 

In this section, the sample derivation of the critical regions for the 
nonlinearity of Figure 7-2 for k ^ 0 and for N = 3 is described. All the 
possible A values A r = (1,0), (1,1), (1,2) are considered. 

The input to the sampler, e(t), is assumed to be a sinusoid as described 

I 

by Eq. 7-17 , and is repeated here, 

e(t) = E cos (W^t + <i>) (7-20) J 

The z -transform of e(t) is 


Ez [ (z - cos 0 ) T) cos <f> - sin w T sin <J)] 

E(z) „ £ £ 

z - 2z cos 0) T + 1 
c 

The period of self-sustained oscillations, T^ is given by 



For N = 3, T = 3T and 
c 


0) 

c 


27T 
3T ’ 


iw T 
~ c 


e j 2 *' 3 


= -0.5 + jO. 866 


( 7 - 21 ) 


(7-22) 


(7-23) 


Then z = e 
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and Eq. (7-21) becomes 

„/ v Ez[(z+.5) cos 6 - .866 sin 61 

E(z) = k | x XL (7-24) 

Z + Z + 1 

i ■■ . 

The waveform for e(t) is shown in Figure 7-5 i f or 6 = 0. Figure 7*6 shows the 
corresponding waveforms of v*(t) for A = (1,0), (1, 1) and (1, 2). 

For A = (1 , 0) ; 

The z-transform of v*(t) is 


V(z) - [M + k(E cos 6 "D ][1 + z" 3 +‘z“ 6 + ...] 


= [M + k(E cos 6 ~D)] <■—-) 

z -1 

The restrictions on E for this mode to exist are 
-30 <6 10° E cos 6 > D 

E cos (60 + (j>)< D 

0 < 6 < 30° E cos 6 > D 

E cos (60 - <j>) < D 

which lead to the following upper and lower bounds on E 


(7-25) 


(7-26) 


-30°< 6 < -0 


D 


E . = . 
mm cos6 


D 


"'max cos (60+6) 


0 < 6 < 30‘ 


min cos 6 


D 


u max cos(60~6) 


(7-27) 


For A = (1,1) : 

The z transform of v*(t) is 


V(z) = [M + k(E cos6 - D) ] [— r — ] 

z -1 


- [M + k (E cos (60 + 6) “ D)] Hi ] (7-28) ‘ 

z -1 
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with the following restrictions on E: 

-30° < <p < 0 E cos (60 + $) > D 

E cos (60 - <j>) < D 
-60° < <j> <-30° E cos <p > D 

E cos (120 + <J>) < D 
The corresponding bounds on E are 


-30° < <p < 0 


E , = D 

min 


cos (60+ $) 


max 


cos (60— <J>) 


-60 < <p < -30 c 


E . 
min 


D 


cos <j> 


E - D 
max 


cos (120+$) 


For A - (1,2) : 

The z- transform of v*(t) is 


V(z) = [M + k (E cos $ - D) ] [— o ] 

z -1 


- [M + k (E cos (60 - <J>) - D) ] [— — ] 

z-1 


- [M + k (E cos (60 +'$).- D)] [ 


3 i 

z -1 


(7-29) 


(7-30) 


(7-31) 


] 
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with the restrictions on E being 

-30° < <t> < 0 E cos (60 - <j>) > D 

0 < <J> <30° E cos (60 + $) > D 

and the corresponding bounds on E being 


(7-32) 


-30° < (j> < 0 


mm 


D 

cos (60— 4>) 


0 < <P < 30' 


E = 00 
max 

E . = D 

mm 


cos (60+<j>) 


E = 00 
max 


(7-33) 


Using E(z) from Eq. (7-24) and V(z) from Eqs. (7-25) ,(7-28) and (7-31) 


the expressions for - 


N(z) 


are as follows : 


For A = (1,0) , 

1 _ _ E(z) _ -E(z-l) ((z+. 5) cos^ - .866 sin 4>) 

N(z) V(z) z 2 fM + k(E cos 4>-D)] 

For A = (1,1) , 

- — ^ = -E(z-l) ((z+,5) cost}) - .866 sin <p) 

Z 2 [Mfk(E cos <p- D)] - [Mfk(E cos (60+4») -D) ] 


For A = (1,2), 


(7-34) ; 

! 


(7-35) 


1 

N(z) 


-E (z-1) ( (z+. 5) cos<|) - .866 sin 4>) 

2 

z [Mfk(E cos <J>-D) ] - z[Mfk (E cos (60-<f>) -D) ] - [Mfk (E cos (60+<J>)-D) ] 


(7-36) 
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With z = -0.5 + j0.866, Eqs .(7-34 ) ,(7-35 ) and (7-36) can be written as 


JL _ 1.5E /-120+<t> 

N (z) ~ (bm + cE) + j(dm + fE) 


where m = M - kD, and the coefficients b, c, d, f are defined in the following 
table. 


Table 7-1 



! 


d 

0.866 

.866 

0 . 

f 

0.866k cos(j> 

0.866k cos<}> 

0.866k (cos<j> + cos (60-<J>)) 


The critical region for each A is the area enclosed by the - ■ / ■ locus 

N (z) 

in the gain-phase plane, as E and <j> are varied over their admissible ranges 
(defined by Eqs. (7-27) ,(7-30) or (7-33) f or that A). 

The procedure used to generate a critical region is as follows: 

1. Choose a value of <J> in the admissible range. 

2. Using this <p in Eq. ( 7- 37) determine the extremums of the 

magnitude and phase of - — (z) as a function of E, with E 

restricted to its admissible range. 

3. Plot these extremums in the gain-phase plane and repeat the 

procedure for all possible 4>’s. The region generated by this 
plot will be the critical <}> region. Figures . 7-14 ,7-15 and 7-16 ' 

show the critical regions for N = 3, A = (1, 0), (1, 1) and (1, 2), 
respectively, for several values of M, k, and D. 
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It should be noted that the extremums necessary in step 2 
above can in fact be easily obtained by use of the upper and lower bounds of 
E [Eqs. (7-27) ,(7-30) and (7-33)]. Also, the extreme magnitude and extreme 
phase occur simultaneously. This result is now proved by use of the expression 
in Eq. (7-37) . 

The magnitude of - as it is defined in Eq. (7-37) ; 


1 


■ N 

l 


r 


2.25E 


rt-cE) 2 + (dm+fE) 2 
Since it is easier to work with 

also the extremums of 

,2 


1 

2 


, the extremums of 


(7-38) 


(which are 


Let F = 


2.25 E 


) will be determined. 
2 


(bm + cE) 2 +(dm + fE) 2 


(7-39) 


If F has any extremums, they must satisfy the following necessary condition: 


3F = 0 = 4.5mE[(b 2 +d 2 )m + (bcKif)E] 

[(bm+cE ) 2 + (dmrffE ) 2 ] 2 


(7-40) 


which gives 
E = 0 


or 


E = - 


2 2 

(b +d )m 
bc+df 


(7-41) 


2 . , 2 , 


In the above results, E = 0 represents a minimum, while E = -(b +d )m/ 

(bc+df) represents a maximum. Let this latter value of E be denoted as E . 

x 


With E = E , the value of F = 
x 

F _ 2.25(b 2 +d 2 ) 
x (cd-bf) 2 


is 


(7-42) 
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For the admissible ranges of <j>, it is seen that the following properties 
of the coefficients are true for all three values of A: 
b > 0 

c > 0 ; 

d >_ 0 

f > 0 (7-43) i 

Using Eq. (7-41) 

E < 0 if m > 0 (M > kD) 
x 

E x > 0 if m < 0 (M < kD) (7-44) 

By physical considerations, the upper and lower bounds of E, which also define 
the admissible ranges of E, are always positive. Thus E^ < 0 is only one of 
academic importance. 

Figures 7-7 and 7-8 show the typical behavior of F, with positive and 
negative values of m, respectively. The asymptotic values of F, as E + + °°, 
are described by 




C7-45) 


Also, it can be shown that F^ >_ F^ for all values of b, c, d, f. 

With N = 3 and any A, it can be shown that E^ is less than the upper and lower 
bounds of E. Thus, for the admissible ranges of 4>, F> and, consequently, 
are always monotonic in E. The maximum and minimum values of - -jjj- are, 
therefore, obtained by use of (i) the allowable maximum and minimum values of E, 
respectively, if m < 0, or (ii) the allowable minimum and maximum values of E, 
respectively, if m > 0. 
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The phase of - ^ from Eq. (7-37);is 


Let g = 


e = 

dm+fE 

bm+cE 


- £ = - 120 ° + 
N 


-1 -dm+fE. 

- tan fes’ 


(7-46) 


-1 


and ip = tan g , then, at constant <f), the maximum and minimum of the 


phase of - — or 0 is determined by the minimum and maximum, respectively, of 
For 4 1 in the range +90°, this is determined by minimum and maximum of g. 
Thus, 0 has only one extremum, and this is when 


E " E y = - 7 1 (7-47) 1 

Figures 7-9 through 7-12show several typical plots of g with m > 0 and m < 0. 
With N = 3 and any A, it can be shown that E is less than the upper and lower 

y 

bounds of E. Thus, for the admissible ranges of <j>, 0 is monotonic in E. The 

1 

maximum and minimum values of the phase of - — are, therefore, also obtained 

by use of the maximum and minimum values of E. The roles of maximum and minimum 

d £ 

can be reversed, depending on the sign of m and the sign of (- ). 

b c 


The above results can be extended to all the cases for N > 3, if the 
corresponding expressions are appropriately modified. The number of modes and the 
ranges of 4> will of course be different for different N. 

The expressions for - l/N(z) for N = 2 through N = 8 have been derived and 
are tabulated in Tables 7- 2 through 1 7_ 8,, 

The critical regions for N = 2 through N = 8 are drawn in Figures 7-13 
through 7-37, for several combinations of k, M, and D. 
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Examination of these critical regions leads to the following conclusions: 

(1) Unlike the critical regions for relay with dead zone (k=0) , 
the critical regions for the nonlinearity in Figure 7-2 cannot 
be normalized in amplitude with respect to D/M. 

(2) The critical regions bounded by the loci of - "• and - 

min 

are symmetrical about the -180-degree axis. The maximum span 
of the horizontal width of these regions is given by 27 rT/T^ = 

2ir/N for N = all even integers >_ 4, and ttT/T^ = tt/N, for N = 
odd integers 3. When N increases the critical regions 
become narrower, and finally approach to a straight line 
along the -180-degree axis when N = “ . 

(3) Since the objective of a great majority of system design 
problems is to arrive at a stable system, for a given system 
whose transfer function G(z) is plotted in the magnitude 
versus phase plot, it is necessary and sufficient that for a 
given sampling period T none of the G(z) points coincides 
with the corresponding critical regions. 

In general, for a quick check on stability, rectangular regions bounded 
by the widths given in item (2) may be used as a necessary condition for 
stability. However, it should be noted that the nonlinear system with k 0 
has less chance of going into self-sustained oscillations than the case when 
k = 0. 

The loci for G(z) for the four cases of the LST system model are plotted 
in Figures 7-38 through 7-47 for N = 2 through 8. The transfer functions for 
these models are tabulated below. In all these cases, the expressions for G (s) 
are given. Then G(s) = G^qCs^^s) is expanded by partial fraction expansion into 



the following form: 



The z-transform of G(s) is written 
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G(s) “ (1 - e Ts ) 


. a _3 + ^2 ■ + \_ 

a 3 2 T s s+b 

s S 


C 1 d l 

+ i^r + iTd 


(7-48) 


where 


G h0 (s) 


1-e 


-Ts 


(7-49) 


G(z) » s[G(s) ] 


(1-z 1 ) z 


f a 3 a 2 a l ^1 C 1 d l 

— j. _ 4. + — £ + 4. — i — 

2 2 s s+b s+c s+d 

s s ) 


(z-1) 


^a.T (z+1) a T a b c 

— 4. — = 4. i. 4. i 4. — 

3 , ,,2 z-1 -bT -cT 


l 2(z-l) (z-1)' 


z-e 


z-e 


z-e 


-dT 


(7-50) 


Case I. 


G (s) 
e 


(Kj s + K q ) 0 ) c 581.5023 (s+4.2) 


Is 2 (s+oj c ) s 2 (s + 70.7) 


Case II. 


G e (s) 


( K l s + K q )K g 82249 (s+4.2) 


9 9 9 9 / 

Is (s +g^s+gQ) s (s +I41.4s+10 ) 


Case III. 


Gg( s ) 


(K n + K u) )s + K A u> 


_0 i_& 


O 


2 w 

Is (s+w ) (s+to ) *" 

C g 


19887. 4(s+3. 6842) 
s 2 (s 2 +10G.7s+2121) 


Case IV. 


G (s) 
e 


(K n + K w )s + K_ to 


-0_ .. 1 . S 




Is 2 (s 2 +g 1 s+g 0 )(s+w g ) G 


(7-51) 


,(7-52) 


(7-53 ) 


281.2928 x 10 4 (s+3.6842) 
s 2 (s 3 +171.4s 2 +14242s+3 X l0 5 ) 



CASE a 3 a 2 a l b l 

I 34.544678 I 7.7363138 -0.10942453 0.10942453 

0.0564239 

II 34.544586 7.7364483 -0.11284781 

-jO. 0017102 

III 34.544632 7.7363682 -0.3835905 0.47625273 

IV 34.544629 ' 7.736475 -0.38701315 0.41178255 


i 


0.0564239 


+ 70.7 


+jO. 0017102 

0 

1 

w_l. 

o 

,» 

*^J 

-0.0926622 

0 

+30 

-0.0123846 

-0.01238466 

+30 


-jO. 0202611 


+j0. 0202611 
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Table 7-2 . 


N = £ 


I _ 


E cos $ 


N(*) 


M + kC B CO$cj>-])') 


c = JB 

L -rnt<n \COS <f>\ 


f = oO 

L *yilAX 


Table 7-3 


N = 3 


N® 


E (Z'0[( Z+OS ^ CoS ^ ~ ° 5l71< ^l 

dZ 2 (w ~kD +fcE cosfy ) + 1>2 ^M-k!D + kE£<?s(^-4>)) 

+ c ( M-kJ)* k E cos (co + fc)) 


1 A 

0- 

b 

c 

i 


E S E 1 

*-»«m j! u ma^ j 

0,°) 

1 

0 

i 

0 

- 3 04 # < O 

$/cOS<j> 

D /cos ( 60 +<p) 

O 30° 

D/ cos ef> 

J) j cos( 60 ~<p) j 

(t,o 

1 

l 

-1 

- 3C% 4> ^ 0 

D/cos(6o+(p % ) 

Djcos (6o -4>) 

_ .. . 

- 60°<cj) ^ -30 

3>/ cos (p 

l)/cOs(l 20 -i-(p ) 

(1,2) 

1 

-1 

-1 

-30°<^ ^ 0 

3> /cos C6o~<f>) 

© o 

0 30° 

J>/ cos (6o+ <P) 

<=o 
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Table 7-4 

N = 4 

I _ ~E ( z. cos<f> - sirufr) 

N(?) a ( M-kD +kE cos 4>) + b (M-fe])-b Esm4>) 


A 

— 

CL 

b 

* 

F 

‘-mm 

P 

t -r»AX 

l\,D 

l 

- 

o 

go 

cos cj> 

3>/-sin <j> 

0$(j>*£4S o 

3)/< CDS <j> 

- - . 

J) J sin (f> 

_ . . 

' 

(2,2) 

1 

l 

-4-S^ (f) ^ o 

V/- Sir » ^ 

cO 



y/ cos cj) 

oO 

. 
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Table 7-5 

N = 5 

| _ - E Lz-')(2 x -22Cosh-Tr/ 5 + \)\_(2-co$iTr/ s )cos<f>-™zv/ s <>iY\Q 

N( 2 ) a (M-JO + If? E cos 4>) ■+ b 2 3 (m- Jo +• k E ccs^^ir/^)) 

+ c z‘ l ^M-l? 3 )+ kE Co.s( 4 > -77-/5-)) 

4* 4 2 ( M-K3>+ kE c + ^l 5 ) 

-1- e ( M-kD + kE c.os(J> - irfe) 



(3,2) 


CO 
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Table 7-6 


N *6 

J _ - £ [ (z -O S)cos<t> - O. 8£6 Str>4>] (z+0 

a'^M + cos4>)z z + b (M -kD + 1 1 ? E cos (Co-t-p^z- 

+ c ( M-k3> +kE cos (&£> -<j>)) 


* 


Vnm 


’mi xx 


-?>0%(bsZO 


O *C (p 30 


y/cos 4> 


3) /cos cj> 


J>/cos(<t>+6o) j 

t 

t 

J)/ CDS L4 > ~&°) 


-30^. <j>^o 


3)/ COs($>t6>o) I n >/ CDS ( (j> -6&) 


- 60^ (p^ - 30° 


- 30 ° 


O ^ 30° 


^>/ COS (p 


3/cos(<P-6o') 


3 )/ CDS (JP+6D) 


oO 


o O 























Table 7-7 


- £ jf ( 2 - CoS-]J-) CX&4> - Sift 4-- $&<?>] (s*- 2 2- CoS^-r + I )(z- l)(l-2z US 4% + 1 ) 
az i> (M-kI> + kECos<t>')+ b z s (M-kD+ l<E ))+ c Z*(M- kD+ !<E u>s(<f> -izy 

+ 4z*(M- kD + kECbs(<t>--^)) 4 - eiCM-kD + ltE cosCd-i-^J) 

-hfi (M- ku + kE cgs(<I> + &)) +9 (m-IcD /<E CosU-Zfj) 



a b c d e f 9 


£3,4) / 



-X 4 (h <0 
14 •’ v 


o* 4 « 


-5^ s0 






■|s*S o 




0 £ ^| 


14 


~ JL < A <-2E 

y M- 


-ft 


0-^-4 


3> 


Cos4> 


r> 


Co s(4> +• ■?£-) 


D 


CO:?c£> 


D 


Cos ( 


cos - 2X) 


D 


cosC<> + ^7) 


X? 


cos(& 


D 


COS (cf> 4 " 7 ^ 





D 


Cos 4 - 23 C) 


T> 


CO s C4 5 4 


D 


cos (. 4 > - 




COS ( 4> ~ S~P) 


x> 


cos C4* 4 


3> 


cos ($4 ^C) 


D 


OS 


p 


COS C4>- ^*) 
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k! = 8 


Table: 7-8 


I __ - £ [(i- 0 . 707 ) CoS <f> -0.7o7Sthfi]U*+j2* + l) 

N A, £m- kD + 1« E cos<£> Jg 5 4 Ce> s (4> + "]£)]& 1 


4- c[M-UD + kEcos ($+£)]$ + 4 £m-M> +ltE.Cos($ -$)J 
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Figure 


7-13 




'HASE (deg) 


7-14 




























PHASE (deg) 


Figure 7-16 















































































PHASE (deg) 


Figure 7-22 




























PHASE (deg) 


Figure 7-24 
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PHASE (deg) 


Figure 7-27 



























PHASE (deg) 

Figure 7-29 














PHASE (deg) 


Figure 7-30 




































PHASE (deg) 


Figure 7-32 















PHASE (deg) 


Figure 7-33 





























PHASE (deg) 


Figure 7-35 
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PHASE (deg 

Figure 7-37 





























PHASE (deg) 


Figure 7-40 


























PHASE (deg) 


Figure 7-42 










PHASE (deg) 


Figure 7-43 






















PHASE (deg) 


Figure 7-45 
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Figure 7-47 
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